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ABSTRACT 


A numerical model has been developed for three-dimensional transient 
conduction based temperature calculations in underwater wet welding on a thick 
rectangular plate. The numerical scheme is based on a fully implicit finite volume 
method. A variable mesh size centered around the moving heat source, and temperature 
dependent thermal properties have been used in the calculations. Convective, radiative 
and boiling surface thermal conditions have also been included. The weld pool region 
itself has been modeled as a solid region of thermal conductivity higher than the 
surrounding unmelted solid region. The validity of the results was checked by comparison 
with Rosenthal's three-dimensional solution for a moving point heat source, and other 
results in the literature. 
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I. INTRODUCTION 


To improve the quality of the underwater welding and to accomplish a reliable, 
permanent underwater wet welding capability has a great importance in today’s industrial 
and military facilities. With the development of underwater wet welding techniques, the 
time and the money required for permanent and temporary repairs of ships and other 
underwater structures can be minimized. Today, the use of hyperbaric welding process 
obtains a limited quality of welding especially for the construction and repairs of 
underwater pipelines. In this process a large pressure chamber is used to keep the water 
away from the workpiece. But, operating this kind of chamber is very expensive and due 
to the limited geometric size, only a few joint configurations can be enclosed in a 
chamber [Ref.l]. The other welding techniques such as double shielding and flux 
shielding also use the water removing theory from the arc area during welding. But, the 
working area must be completely prevented from water for satisfactory welding results 
[Ref. 2]. 

Currently, underwater wet welding is used for the temporary repair needs. 
Because of their poor quality compared to surface (air) welds (they obtain 80% of the 
tensile strength and 50% ductility of the surface welds [Ref. 3].), they must be replaced as 
soon as possible. Therefore, the development of a more efficient wet welding technique is 
the only solution to this problem. 

The surrounding water environment during wet welding causes rapid cooling and 
steep temperature gradients in the weld area behind the arc [Ref.4]. Because of the 
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extremely complex nature of the heat transfer phenomena between heated surface and the 
water environment, a numerical model simulation is necessary. 

In the present study, a numerical model has been developed for transient, three- 
dimensional conduction heat transfer in underwater welding process on a thick 
rectangular plate. The numerical scheme was based on fully implicit finite volume model, 
including convection, radiation and boiling surface thermal boundary conditions. The 
different regimes of boiling were accounted for on the surface. A variable mesh size 
centered around an arc source moving at constant speed was used to determine 
temperature variations inside and around the weld pool. The weld pool region itself has 
been modeled as a solid region of thermal conductivity higher than the surrounding 
unmelted solid region. The input data from the previous studies based on different 
methods were used to check the accuracy and the validity of the numerical method. 
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II. BACKGROUND 


A. PREVIOUS STUDIES 

Rosenthal did the most important early work on the theory of the effect of moving 
sources of heat in the late of 1930s. He studied the fundamentals of this theory and 
derived equations for two-dimensional and three-dimensional heat conduction in a solid 
when a moving source is in use [Ref. 5,6], The analytical solution derived by Rosenthal 
was based on the principle of a quasi-stationary thermal state. The quasi-stationary 
thermal state represents a steady thermal response of the weldment with respect to the 
moving coordinates. In other words, the origin moves with the heat source, and for an 
observer at origin, the temperature distribution and the pool geometry do not change with 
time [Ref. 6,7]. Rosenthal's solution for three-dimensional heat flow during welding is 
as follows [Ref. 6]: 


2 

Q 


= exp 


f U(R-x) ^ 
2a, 


( 2 . 1 ) 


where 

( 2 2 2V /2 

x +y +z ) 

T = temperature, 

T 0 = workpiece temperature before welding, 
k s = thermal conductivity of solid, 

Q = heat input to the workpiece, 
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U - welding speed, 

x = the distance that the heat source traveled, 

a s = thermal diffusivity of solid (i.e., k s / pC s , where p and C s are density and 
specific heat of solid, respectively). 


The assumptions used by Rosenthal to derive the equation above are as follows 
[Ref. 7]: 

1. Point heat source. 

2. No melting and negligible heat of fusion. 

3. Constant thermal properties. 

4. No heat loss from the workpiece surface. 

5. Infinitely wide workpiece. 

Although Rosenthal's assumptions help to simplify the mathematical analysis 
involved, there are some significant deviations between theoretical and experimental 
results. For instance, as a result of the point heat source assumption, the temperature at 
the weld centerline goes to infinity even though the power of the heat source is finite. 
Also, the values of thermal properties change with the temperature and neglecting the 
heat fusion gives considerable errors. [Ref. 7] 

Tanaka did the first studies on the application of the mathematical analysis of 
non-stationary heat flow to the practical problems. Naka and Masubuchi studied the 
mathematical analysis of non-stationary heat flow. They used dimensional expressions to 
make the numerical calculations simple. Nippes and Savage studied the cooling rates of 
heat affected zones by using a graphical approach method. Suzuki found an analytical- 
empirical method of studying the effects of welding parameters and determined the 
cooling rate from welding conditions with the help of a monograph. [Ref.6] 
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Adams derived new equations from Rosenthal’s equation by using the fusion line 
as a boundary condition and calculated the peak temperature at a distance from the fusion 
boundary at the weldment surface [Ref. 7]. His equation for three-dimensional heat flow 
is as follows: 


1 _ 5 . 44 ^ 0 ^ 

T p -T 0 ~ QU 

where all the terms are defined in equation (2.1) except 

T p = peak temperature at a distance Y from the fusion boundary, 

T m = melting temperature, 

Christensen also used analytical-empirical approach and derived dimensionless 
equations based on Rosenthal's three-dimensional equation. He explained the relation¬ 
ships between the welding conditions and the weld bead geometry. [Ref. 7] 

Recently, numerical analysis methods and computer programs have been 
commonly used to develop the previous assumptions. In 1965, The Battelle Institute 
Geneva Laboratory in Switzerland conducted a computer-aided study about analyzing of 
heat flow in weldments [Ref. 6]. The University of Wisconsin and McDonnell Douglas 
Aircraft Company also conducted similar studies in the same year [Ref. 6]. In 1970, 
M.I.T. researchers studied heat flow during underwater welding. They also developed 
finite element programs on heat flow during welding [Ref. 6]. Oreper and Szekeley 
examined stationary, axisymmetric TIG (tungsten-inert-gas) welding process with a 


2 + 


f UY ^ 2 


\ 2a s J 


+ - 


T m ~T 0 


( 2 . 2 ) 
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moving boundary by using finite difference method. Their formulation contained the 
affects of transient conduction, electromagnetic, buoyancy and surface tension forces 
[Ref. 8]. Kou and Wang performed computational studies of the GTA welding process. 
They presented a computer simulation of three-dimensional convection for an arc source 
moving at constant speed. They considered electromagnetic, buoyancy and surface 
tension forces on the pool surface. They found very good agreement between the 
calculated and observed fusion boundaries [Ref. 9]. They also studied the computer 
simulation of three-dimensional convection in laser melting by considering the buoyancy 
force and the surface tension gradient at the weld pool surface [Ref. 10]. Correa and 
Sundell studied axisymmetric stationary arc source by using different grid sizes for 
computation of flow and electromagnetic fields [Ref. 11]. Saedi and Unkel developed a 
thermal-fluid model of the weld pool. Their model was based on the stationary arc. To 
describe the weld pool geometry, they matched the convective and the conductive heat 
fluxes at the weld boundary by using an iterative calculation method [Ref. 12]. Zacharia 
et al. made three-dimensional calculations on the effects of the heat source in the 
stationary GTAW process. He indicated that a depressed area formed at the weld pool 
center because of an outward fluid flow caused by surface tension force [Ref. 13]. Kim 
and Na developed a model on heat and mass flow for stationary, GTAW process with 
electromagnetic, buoyancy and surface tension forces. They used numerical mapping 
method for calculations [Ref. 14]. Ramanan and Korpela compared the effects of thermo¬ 
capillary and Lorentz forces on the flow pattern in a stationary weld pool with buoyancy 
forces. They used multi-grid methods and a local grid refinement technique with 
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axisymmetric stationary arc source [Ref. 14]. Ule, Joshi and Sedy determined three- 
dimensional transient temperature variations in the GTAW process by using the explicit 
finite difference method. They used different mesh sizes and temperature dependent 
thermal properties. They also considered convective and radiative surface thermal 
conditions during calculations [Ref. 16]. Kanouff and Greif studied the unsteady 
development of an axisymmetric arc weld pool in GTAW process. They used moving 
grids to follow the phase change boundary and considered the effects of Marangoni, 
Lorentz and buoyancy forces in the calculations [Ref. 17]. Joshi, Dutta, Schupp and 
Espinosa developed a three-dimensional numerical model to describe the flow circulation 
phenomena in al uminum weld pools under non-axisymmetric Lorentz force field [Ref. 
18,19]. 


B. FINITE VOLUME METHOD 

The finite volume method is one of the simple and well-established 
Computational Fluid Dynamics (CFD) techniques that were originally developed as a 
special finite difference method [Ref. 19]. The stages of the numerical algorithm in this 
method are as follows [Ref. 19]: 

1. Formal integration of the governing equations of fluid flow over all the finite 
control volumes of the solution domain. 

2. Discretisation involves the substitution of a variety of finite difference type 
approximations for the terms in the integrated equation representing flow 
process such as convection, diffusion and sources. This converts the integral 
equations into a system of algebraic equations. 

3. Solution of the algebraic equations by an iterative method. 
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In the control volume integration, the calculation domain is divided into discrete 
control volumes. There is only one control volume surrounding each grid point and the 
boundaries of the control volumes are positioned side to side in the middle of the distance 
between the grid points. Integrating the governing differential equation over each control 
volume derives the resulting discretised equations. The discretised equations express the 
conservation of quantities such as mass, momentum and energy for each control volume. 
This characteristic exists for any number of grid points. The numerical solution insures 
the validity of the conservation principle over the whole calculation domain for the 
related quantities. Versteeg and Malalasekera [Ref. 20] wrote " This clear relationship 
between the numerical algorithm and the underlying physical conservation principle 
forms makes finite volume method much simpler to understand by engineers than finite 
element and spectral methods". To understand the finite volume method better an 
illustrative example can be given for one-dimensional steady state heat conduction 
situation. [Ref. 20,21,22] 

1. One dimensional steady state heat conduction 

The governing equation for one-dimensional steady state heat conduction is 


cbc v dx ) 


+ S = 0 


where 

k = thermal conductivity, 

T= temperature, 

S = the rate of heat generation per unit volume. 


(2.3) 
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The domain is divided into small and nonoverlapping control volumes. A part of 
one-dimensional grids generated is shown in Figure 2.1. Here, The grid point under 
construction is denoted as P and the neighboring nodes to the east and west are denoted as 
E and W respectively. The lower case letters e and w denotes the east and west faces of 
the control volume. The distance between the nodes W and P is denoted by Sx lvp and 

between P and E is denoted by dx PE . The distances between w and P and between P and e 
are given by Sx wP and dx Pe respectively. The control volume width is shown as 
Ax = <fc we 



Figure 2.1 One-dimensional grid [Ref. 20] 


Integrating equation (2.3) over the control volume gives. 


(i dT ) 

k — 
dxJe 



(2.4) 


9 




By assuming a piecewise-linear profile assumption (Figure 2.2), equation (2.4) 


can be evaluated as 


t e -t p 


Sc 


-k. 


f T —T ^ 
1 P L W 


PE J 


Sx 


+ 5Ax = 0 


WP J 


(2.5) 


where S is the average value of S over the control volume. 



Figure 2.2 Piecewise -linear Profile [Ref. 22] 


By arranging the terms, the final discretised equation can be written as 

QpTp = a E T E + o w T w + b 


where 


a 


E ~ 


Sc p E 


( 2 . 6 ) 


(2.7a) 
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(2.7b) 


a P =a E + a w 


b = S Ax 


(2.7c) 


(2.7d) 
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IH. MODEL DEVELOPMENT 


A. DEFINING THE WORKPIECE 

The workpiece has been defined as a thick rectangular plate (Figure 3.1). 



Figure 3.1 The workpiece 


where 

1. Right Lateral Face (East) 

2. Left Lateral Face (West) 

3. Back Face (North) 

4. Front Face (South) 

5. Top Face (Top) 

6. Bottom Face (Bottom) 

Convective heat transfer coefficients for each face have been defined as follows: 

Right-Lateral Face: h e 
Left-Lateral Face : h w 
Back Face : h n 

Front Face : h s 

Top Face : h, 

Bottom Face : h b 

The same method has been also used to define thermal conductivity and heat flux terms 
for each face. 
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B. BOUNDARY CONDITIONS 


1. Top Face 

By adding radiative and convective heat loss from the top face, boundary 
condition equation can be defined as 

~ k %‘ h(T "“ - r »> + «" +<K fc -C) (3.1) 

In equation (3.1), for the radiation term, T mr can be neglected, (j^,, » C).The 

distance between the node point P and the face can be taken as Ax/2, where the distance 
between the node point P and the neighboring nodes is Ax . Again, for the radiation term, 
because of the small distance between the node point P and the face, by assuming 
T wotl = T P , we have 


-k 


(TwaU Tp ) 

Ax/2 


= KT waII -TJ + q"+asTp 


(3.2) 


Equation (3.2) can be opened as 


2k 

Ax 


[ wall 



= hT wa ,i -hT K + q"+ crsTp 


(3.3) 


Equation (3.3) can be arranged as 
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(3.4) 


L wall 


(2k ' 

— + h 
\&x , 


~T P + hT x -q"+ <jsTp 
Ax 


Equation can be written as 


[ wall 


Ik 

—Tp +hT n -q"+(jsTp 
Ax 


2k 
— + 
Ax 


h 


(3.5) 


The unknown temperature for the top face points is found as 


-wall 


2kT P + hAxT a —q" Ax - crsAxT* 
2k + hAx 


(3.6) 


where 


<l",op = <l" sourced" ^'boiling 


(3.7) 


q"source = arc source heat flux, 

q'\ = a constant arbitrary heat flux which may be applied to the top face, 

q"tom ng = foiling heat flux, 


2. Other Faces 

By using the same method from the top face (without radiation). 
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(3.8) 



2kTp +hAxT , -g"Ax 
2A: + AAx 


3. Simulating the Arc 

To simulate the heat input from the arc to the workpiece, it is assumed that the 
heat input distribution of the arc have a Gaussian distribution on the top face of the 
workpiece. The general equation is 

00 

Q = q 0 Je r ° 2nrdr (3.9) 

0 

where 


Q - the total heat input into the workpiece, 
q 0 = the volumetric energy generation rate, 

r 0 = the radius of the heat input distribution, 
d = the exponential factor, 

By solving equation (3.9), the volumetric energy generation rate can be found as 


#0 ~ 


Qd 


7U\ 


(3.10) 


4. Boiling Heat Transfer 

Modes or regimes of boiling and the related equations can be classified as follows 
(where AT e =T S - T sat ): 
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a. 


Free Convection Regime (A T e <5 °C) 


In this regime, natural convection effects determine the heat transfer 
between the heating surface and surrounding liquid. Recommended correlations for upper 
surface of heated plate are [Ref. 25] 


Nu l = 0.54 Ra" 4 

(lO 4 <Ra L <10 7 ) 

(3.11) 

Nu l =0.15 Ra ] / 3 

(lO 7 <Ra L <10") 

(3.12) 


where the Rayleigh number, 


R ZW.-T.V 

va 


(3.13) 


here 


g = gravitational acceleration, m/s 
l3 = volumetric thermal expansion coefficient, K' 1 
v = kinematic viscosity, m /s 
a = thermal diffusivity, m 2 /s 

L = characteristic length, L = Plate surface area (A s )/Perimeter (F) 


7 Ntnk 

h = - 

L 


(3.14) 


and the value of heat flux is, 
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q"= h(T s -T x ) 


( 3 - 15 ) 


b. Nucleate Boiling Regime (5 °C < A T c < 30 °C ) 

The most useful nucleate pool boiling correlating equation was developed 
by Rohsenow [Ref. 23, 24, 25,], 




g(Pi -pS 

1/2 

c,AT, ) 

cr 


Pr " ) 


(3.16) 


where the subscripts s, 1, and v express surface, saturated liquid state and vapor state. The 
definition of each term in equation (3.16) are as follows: 


fj. t = viscosity of the liquid, kg/ms 

h fg = latent heat of vaporization, J/kg 

g = gravitational acceleration, m/s 2 

P' = density of the saturated liquid, kg/m 3 

P v = density of the saturated vapor, kg/m 3 

° = surface tension of the liquid-to-vapor interface, N/m 

c pl = specific heat of saturated liquid, J/kgK 

A T e = T s -T sal 

Pr, = Prandtl number of the saturated liquid 
n =1.0 for water, 1.7 for other fluids 

C s J = empirical constant that depends on the nature of the heating 

surface fluid combination and whose numerical value varies 
from system to system 


But, in underwater welding, the surrounding water temperature is below the 
saturation temperature (between 0 °C and 30 °C ). This is called as the heat transfer to a 
subcooled liquid. For the subcooled boiling, the heat flux can be estimated as [Ref. 23] 
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0 = 0 . 


(3.17) 


r 

< 1 4- 

2^/ (T sat T liquid ) 

24 

r 2 i 

Pv 

1/4 " 

i i 

■yjna t z 

nh fgPv 

_og(p,-p v )_ 



where 

k 

1 = thermal conductivity of the liquid, W/m.K 
a ‘ = thermal diffusivity of the liquid (k / pc p ), m 2 /s and. 


t = 



G 

1/2 

" 2 1 
Pv 

_g(j>l -Pv). 


_og(p ! -p v )_ 


(3.18) 


c. Transition Boiling Regime (30 °C < A T e < 120 °C) 

For the transition-boiling regime, no sufficient theory has been derived. 
This regime is between the maximum and minimum heat fluxes where [Ref. 23,25], 


0"max = O-149 h fg p x 


°g{Pl ~ Pv) 


nl/4 


(3.19) 


0"min = 0-09 p v h fg 


gajp, -p v ) 
. (pi+PvY 


—11 / 4 


(3.20) 


By assuming a linear heat flux distribution in the transition boiling regime 
(Figure 3.2), it can be written, 


lOg(0 M max )-lOg(0"min ) __ l0g(g") ~ lOgCg"^ ) 

Iog30-logl20 log AT e -log 120 


(3.21) 
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by arranging equation (3.21), 


log 


log(<7") = • 


( Q" \ 

“ max 

J 

( 1A \ 


AT 


log 


30 

Vl20y 


120 


log T-r +log(^" min ) 


(3.22) 


and the heat flux at a point in transition boiling regime can be found as 


q"= io log(? ") 


(3.23) 



Figure 3.2 Linear Distribution of Heat Flux in Transition Regime 

d. Film Boiling Regime (AT e > 120°C) 

In film boiling regime for flat horizontal surfaces, Westwater and Breen 
recommended the following correlation for conduction heat transfer coefficient [Ref.23], 
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he = 0.59- 


g(p, -p;)pX 3 [h fg + 0 - 68 c^Ar.]] 


^vAr e 


1/4 


where 

^ = viscosity of the vapor, kg/ms 

K = thermal conductivity of vapor, W/mK 

c pv = specific heat of saturated vapor, J/kgKand, 


X = 2m 


a 


, 1/2 


g(Pl-Pv)] 


(3-24) 


(3-25) 


Bromley suggested combining conduction and radiation heat transfer coefficients 
[Ref. 23], 


htotal — hc+0.75hr 


where 



f r r A > _ \ 

1 s 1 sat 

T s ~T sat , 


here 

* 

s s = surface emissivity 

T s = absolute surface temperature 

and the resulting heat flux in this regime is found as 


q"~ hwtai A T e 


(3:26) 


(3.27) 


(3.28) 
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C. COEFFICIENTS USED IN THE EQUATIONS 


The coefficients used deriving the discretised equations are as follows: 


k„A„ 




Sx 


PE 


a w ~ 




Sx 


ap 


a N = 


kA 

n ti 


5y 


PN 


M 

fysp 


a = - J 
u s 


a-r = 


M, 


r & 


PT 


a B ~ 


&BP 


o AV 

a p = pc - 

P At 


a p = ®{a E + a w + a N +a s +a T + a B )+ a 


coeffe = *'** 
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fliacoeffn ■■ 


fypN 


2k n + h n dy PN 


fluxcoefft = 


dz 


PT 


2k, + h,dz PT 


radcoeff = as ^PT 


Ik, + h,dz PT 


fluxcoeffs = 


fysp 


2k s + h s 5y sp 


flvxcoeffb = 


dz 


BP 


2k b + h b dz BP 


D. DERIVATION OF THE EQUATIONS 

The problem is governed by the equation 



dx J 



dT* 
dy j 


d (. dT\ 
+ — k — 
dz v dz J 


1. Interior Nodes 


(3.29) 
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Integration of equation (3.29) over the control volume and over a time interval from t to 
t + At gives 



(3.30) 


This may be written as 




Adxdt + 


t+At 


!!- *- 


Adydt 


t+At 




/ cv 


ar 

dz ) 


Adzdt 


(3.31) 


where A is the face area of the control volume and dx, dy, dz are the dimensions of the 
control volume. Equation (3.31) may be written as 
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(3.32) 
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By assuming the grid point value of the temperature at a node to prevail over the whole 
control volume, the left hand side of equation (3.32) can be written as 



(3.33) 


where T p refers to temperatures at time t and T p refers to temperatures at time t + At. 
By applying central differencing to the diffusion terms on the right hand side equation 
(3.32) can be written as 
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(3.34) 


The values of T p ,T e ,T w ,T n ,T s ,T t and T B vary with time. The time integral can be 

calculated by using temperatures at time t or at time t + At or, a combination of 
temperatures at time t and t + At. This approach may be generalized by defining a 
weighting parameter © between 0 and 1. 



t+At 

ft 


dt = \-jT r + (l - Bft Jl/ 


t 


(3.35) 
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In equation (3.35) for 0 =0 the temperature at time level t is used (Explicit scheme); for 
0=0.5 the temperatures at t and t + At are equally weighted (Crank-Nicholson scheme); 
for © =1 the temperature at time level t + At is used (Fully implicit scheme). By using 
equation (3.35) and dividing At throughout, equation (3.34) may be arranged as 
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(3.36) 


Equation (3.36) may be re-arranged as 
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(3.37) 
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By putting the defined coefficients from the previous section, equation (3.37) may be 
written as 

[a® + ®{a E +a w +a N +a s +a T +a B )]r p = a E [®T E + (l - ®)T e ] 

+ a w [®T W + (1 - ©W ]+ a N [®T n + (1 - ©)T W 0 ] + a s [©r s + (l - ©fc 0 ] 

+ a T [©7V +(l - ©)r r 0 ]+ a B [©r B + (l - ©)r 5 0 ] 

+ [a® - (1 - ©)(a £ + a w + a N + a s + a T + a B )} T P (3.38) 


Finally, by grouping the known and the unknown terms at each side, the discretised 
equation is found as 


a p T P -®a E T E -®a w T w -®a N T N -®a s T s -@a T T T -Qa B T B = (1-0) 

[a E Tl +a w T* + aX + a s T s + + ***?]+ l a l ~( l ~ 0 )<>£ + 

+ a N + a s +a T + a B )\r p (3.39) 


2. Left-Front-Top Corner 
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By using the same method from equations (3.30-3.33), equation (3.29) can be written as 
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(3.40) 


W 


In equation (3.40), T w and T s take the value of equation (3.8) and T t takes the value of 
equation (3.6). By using these values and equation (3.35), we have 
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(3.41) 


Equation (3.41) may be arranged as 
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2kA (hsM* zl± 2 M \ h &lil 
3y S p 1 + h s dy SP J & PT 


-q'\ op &pt ~crs& PT Tp 
2k, t +h t & PT 


(3-42) 


Equation (3.42) may be re-arranged as 


h A«- juf Mk s 

A/ + h w &c WP J Sy PN fy$ P K 2k s +h s 5y SP j 
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K&ptT^ _ ? top &pt a£&p T T F ^ | h v 8x iV i>T m 

2k t + h,& PT 2k. +h t & PT 2k, + h,Sz PT j v 2A' u + h w 5x lVj: 

V\ K? ) , 2kA ( Kfyffff )| 2 ^ 4 f 

2^k J fysp k 2^j + V*V 2& f +h s $ysp > ^pt v 2& ( + h,& PP 


‘ftop fepT ae5z PT T p \ 
2k, + h,Sz pr 2k, + h,Sz PT 


By putting the defined coefficients from the previous section, equation (3.43) may be 
arranged as 


[ a p + ©[a zr + 2 <% {coeffw) + a N + 2 a s {coeffs) + 2 a T {coefft) + a B ])r p 
= a £ [&T e + (1 - Q)T° ]+ a N [&T n + (1 - ©)r° ] + a, [er fi + (1 - ©)r s ° ] 

+ [ a °p -0- ©)[«£ + 2<% {coeffw) + a A , + 2u s (coe#s) 

+ 2a r (coefft) + a B ]fr p + ©[2%. [{coeffwff^ - {fluxcoeffw)q" J 
+ 2a s [{coeffs)ff - {fluxcoeffs)q" s ] 

+ 2a T \coefft)T x - {fluxcoefft)q'\ op -{radcoefft)T p ] 

+ (1 - ©)[2a„. [{coeffw)T w - {fluxcoeffw)q\] 

+ 2 a s {coeffs)T K - {fluxcoeffs)q" s ] 

+ 2 a T \coefft)T x - (fluxcoefft)q” top ~{radcoefft)T P ] (3.44) 


by putting heat flux terms from equation (3.7), equation (3.42) may be written as 

[a° P + @[a E + 2 a w {coeffw) + a N + 2 a s ( coeffs ) + 2 a T ( coefft ) + a B Jr, 

= a E [er £ + (i - ®)T° ]+ [©r, + o - ©)r„° ]+ a, [er, + (i - ejr," ] 

+ [a? - (1 - ©)[<>£ + 2 a w ( coeffw ) + a,. + 2a. (coeffs) 

+ 2 a T {coefft) + a B Jr/ + [2a w {coeffw) +2a s {coeffs) + 2 a T {coefft )Jr oo 
- 2a w {fluxcoeffw)q" w -2a s {fluxcoeffs )q" s -2 a T {fluxcoefft )q" t 
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- 2 a T {fluxcoeffi )q" source -2 a T (fluxcoefft)q" b0lling -2 a T {radcoeff)Tp 


(3.45) 


and, the discretised equation may be written as 


[a° P + ®[a E + 2 a w ( coeffw ) +a N + 2 a s {coeffs ) + la T (coefft)+ a B IJZ> 

-@a B T B -&a N T N -&a,T t =(l-©)[a j r i » +a„T’ +a B T°] 

+ a° p - (1 - ©)[a £ + 2 a w {coeffw) + a N + 2 a s {coeffs) + 2 a T {coefft)+ a B Jr p ° 

+ [2a w {coeffw) + 2 a s {coeffs)+ 2 a T (coefft)]T m - 2a w {fluxcoeffw)q'\, 

- 2a s {fluxcoeffs)q" s -2a T {fluxcoefft)^',-2a T {fluxcoefft)q” source 

- 2 a T {fluxcoefft)q" boiling -2 a T {radcoeff)T A p (3.46) 


3. Front-Left Edge 



Figure 3.5 Front-Left Edge 


By using equations (3.28-3.31), equation (3.27) can be written as 
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(3.47) 


In equation (3.47), T w and T s take the value of equation (3.8). By using these values 
and equation (3.35), we have 


pc— — -^-AV =0| 


At 


k - A -(T c -T F ) 2KA ~ 


Sc 


PE 


Sc, 


WP 


rp JKTp +KSCmT x cf w Sc wp 

J p 




2k w + h w Sc K7> 


J 


k„A 2 k'A, 

+i fhL (T N ~T P )—^ 

fypN WSP 


S S \ T _ 
L p 


KA 


BP 


(T P -T b ) 


+(!-©) 


k„A„ 


2 K T p +h s fy sp T„ -q" s fy SP 
2 k s + h s Asp 
( 


+| ~(T t -T r ) 


PT 


Sc 


/j>0 j 0\ f j-0 ^vTp jAAgA w 

' E P* e. P . 7 cv 






w v 


24 +^ H < 5 C([/p y 


k„A/jo n<>\ 2Jc s A s 2k s T% +h s fy SP T x q s fy SP 
+ c v / at - t W c ^/> r: rr 

9W 0V l, 2k s +h s fy SP 


k,A, 


k k A k 


dz 


(T°-Tph-^-CTp -Tjj) 


PT 


dz 


BP 


(3.48) 


Equation (3.48) may be arranged as 
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Equation (3.49) may be re-arranged as 
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(3.50) 
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By putting the coefficients from the previous section, equation (3.50) may be written as 


[a° P + e[a E + 2 a w (coeffw) +a N + 2 a s (coeffs) + a T +a B Jr, 

= °e¥ t b +( 1 -@)T° +fl*[©r„ +(l-e)7-„“]+a r [® 7 ’r +(l-®)7?] 

+ a B [ 07 J + (1 - ]+[aj -(1 - ©)[a 5 + 2a w (coeffw) -a,, 

+ 2 a s (coeffs) + a T + a B J7>° + ®[la w {(coeffw)ff - (fluxcoeffi)q" w ] 

+ 2 a s \(coeffs)T x - (fluxcoeffs)q n J+ (1 - ©)[2 a w [(coeffw)T K - (fluxcoeffw)q" w ] 

+ 2a s [(coeffs)ff - (fluxcoeffs)q" s ] (3.51) 


Equation (3.51) may be arranged as 


[a£ + ©[fl £ + 2 a w (coeffw) +a N + 2 a s (coeffs) + a T +a B 
= a E [er £ + (1 - ©)T° ]+ a N [qt n + (1 - ©)T° ] + a T [©r r + (1 - ©)T t ° ] 

+ a B [©^fi + (1 - ®)Tb - (1 - ©)[«£ + 2 a w (coeffw) +a N 
+ 2a s (coeffs) + a T + a B Jr/ + [2 a w (coeffw) + 2 a s (,coeffs )]T 00 
- 2 a w (fluxcoeffw)q\ -2a s (fluxcoeffs)q ", (3.52) 


Finally, the discretised equation may be written as 


[a° r +&[a c + 2 a w (coeffw) +a„+ 2a s (coeffs) + a T + a, Jr, - ®o £ r £ - ®a„T„ 

- @a r T r -%a„T B = (1 - @)[a s T° + + °X]+ [°° r -0 -©)k 

+ 2 a w (coeffw) + a N + 2 a s (coeffs) + a T + a B Jr/ + [2a ;v (coeffw) + 2 a s (coeffs)^ 
- 2a w (fluxcoeffw)q" w -2a s (fluxcoeffs)q" s (3.53) 




4. 


Bottom Face 
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Figure 3.6 Bottom Face 


By using equation (3.30-3.33), equation (3.29) may be written as 
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(3.54) 


In equation (3.54), T B takes the value of equation (3.8). By applying this and equation 
(3.35) to equation (3.54), we have 
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(3.55) 


Equation (3.55) may be arranged as 
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(3.56) 


Equation (3.56) may be re-arranged as 
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=-^%r £ +(i-0)r”]+%t[ err +( i_0 ) r;]+M-[®7-„ +(i-©)r„°] 

COCpE OX ]Vp fypN 

+07V + a - e)r s ° ]+M- [®r r +(i - ©)r r ° ]+ - a - ®{-|^ 

^sp &pj- A/ |_ atj^ 

! | Kt^n | | | ^b^b K^jBP j'O j 0 2k b A b 

fypN 3ySP dzPT &BP Kp'^b + BP J _ BP 

K^BP^x _ Q b ^BP + ( 2 _ 0 ) ^Jh^b hb^BpTta 

\2k b + h b dz BP 2k b + h b Sz BP j &bp \2k b + k b Sz BP 


tf'b &BP 

2k b + h b Sz BP J_ 


(3.57) 


By using the defined coefficients, equation (3.57) may be written as 


[a° P + ®[a E + a w +a N +a s +a T + 2 a B (coeffb)^r p = a E [©T £ + (1 - ©)r £ ° ] 

+ a w [®T w +(l-©)^]+a„[©r„ +(l-®)7^ 0 ]+fl s [©r 5 + (1-0)T S °] 

+ a r [©r r + (1 - &)T °]+ [a° P - (1 - ©)[a £ + a w +a N +a s +a T + 2a B (coeffb)^r p 
+ ® [2a B [{(coejjfb)T x ) - ({fluxcoeffb)q " b )] 

+ (1 - ®)\la B ]^coejfb)T x ) - {(fluxcoeffb)q" h )]] (3.58) 


And, the discretised equation may be written as 


[a° p + &[a E + a w +a N +a s +a T + 2a B {coeffb)^ p - ©a E T E - © a w T w - © a N T N 

- © a s T s - © a T T T = (1 - ©)[a E T E + a w T° + a N T° + a s T° + a T T P ]+ [a p 
~ (1 - ®)[a E + a w +a N +a s +a T + 2 a B (coeffbj§T P + [2a B (coeffb)]T x 

- 2 a B (fluxcoeffb)q" b (3.59) 


The resulting discretised equations for the other parts of the workpiece are as follows 
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5. Left-Back-Top Corner 


[a° p + ®[a E + 2 a w ( coeffw ) + 2 a N (coeffri) + a s + 2a r ( coefft) + a B Jr, -®a E T E 

- ® a s T S - Qa B T B = (1 - ®i?E T E + a S T S + a B T B ]+ [«£ “ 0 ~ <=>)[«£ + 2 % (. COefjfw ) 
+ 2 a N (coeffri) + a s + 2 a T (coejft) + a B Jt p ° + [2 a w (coeffw) +2 ct N (coeffri) 

+ 2a T (coefft)]T x - 2a w (fluxcoeffw)q'\ -2a N (fluxceeffii)q" n -2a T (fluxcoefft)q", 

- 2 a T (fluxcoeff)q" source -2a T (flvxcoeffi)q" boiling 

- 2 a T (radcoeff)Tp (3.60) 


6. Right-Back-Top Corner 

[a° p + ®\2a E (coeffe) + a w + 2a N (coeffri) + a s + 2a T (coefft) + aJfT p - ®a w T w 
-®a s T s -®a B T B =(l-®)[a w T°+a s T s ° +a B Tf]+[a° p -(l-®)[2a E (coeffe) 

+ a w + 2a N (coeffri) + a s + 2 a T (coefft) + a B ]T P + [2a E (coeffe) +2 a N (coeffri) 

+ 2 a T (coefft)]T x - 2a E (fluxcoeffe)q" e -2a N (fluxcoeffh)q '\ n -2 a T (fluxcoefft)q"\ 

- 2 a T (fluxcoefft )q" source -2 a T (fluxcoefft)q 

boiling 

- 2 a T (radcoeff)Tf (3.61) 


7. Right-Front-Top Corner 

[a° p +®[2a E (coeffe) + a w + a N + 2 a s (coeffs) + 2 a T (coefft) + a B f P -®a w T w 
~ ®a N T N - ®a B T B = (1 - ®)[a w T„ +a N T° + a B T B ]+ [a° F -(1 - ©)[2 a E (coeffe) 
+ a w + a N + 2 a s (coeffs) + 2 a T (coefft) + a B ]r p + [2a E (coeffe) +2 a s (coeffs) 

2a T (coefft)^ - 2a E (fluxcoeffe)q'\ -2a s (fluxcoeffs)q n s -2a T (fluxcoefft)q'\ 

- 2 a T (fluxcoeff)q" sourc -2a T (fluxcoefft)q 

boiling 

- 2 a T (radcoeff)T p 


(3.62) 
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Left-Back-Bottom Corner 


[a° p + ©[a £ + 2 a w {coeffw ) + 2a N {coeffh) + a s +a T + 2 a B (coejfb)§r p -®a E T E 
~ ® a s T s ~®a T T T = (1 - ®)[a E T E + a s T B + a T T E \+ [a° p -(1 -©)[a £ + 2 a w {coeffw) 
+ 2 a N (coeffh) + ct s +a T + 2a B {coeffbj\Tp 

+ [2 a w {coeffw) +2a N {coeffh) + 2 a B {coeffb)]T x - 2a w {fluxcoeffw)q" w 

-2 a N {flwccoeffh )q" n -2a B {flvxcoeffb )q'\ (3.63) 


9. Right-Back-Bottom Corner 

[a° p + ®\2a E {coeffe) + a w + 2a N {coeffh) + a s +a T + 2a B {coeffbj§T P - ®a w T w 
-Qa s T s -®a T T T = (1 - ®)[a w T° + a s T° + a T T T °]+ [a° p - {l-®)[2a E {coeffe) 

+ a w + 2a N {coeffh) + a s +a T + 2 a B {coeffb)]Tp 

+ [2a E {coeffe) +2a N {coeffh)+ 2a B {coeffb)^ -2a E {fluxcoeffe)q" e 

- 2 a N {fluxcoeffh)q" „ -2a B {flvxcoeffb)q" b (3.64) 


10. Left-Front-Bottom Corner 

[a° p + ®[a E + 2a w {coeffw) + a N + 2a s {coeffs) + a T + 2a B {coeffb)^ff p - © a E T E 
-®a N T N -®a T T T = (1 - ®)[a E T E +a N T° + a T T°]+[a° p -(1-©)[a £ + 2 a w {coeffw) 
+ a N + 2a s {coeffs) + a T + 2 a B {coeffb)]Tp 

+ [2 a w {coeffw) +2a s {coeffs) + 2a B {coeffb)]T x - 2a w {fluxcoeffw)q'\ 

- 2 a s {jluxcoeffs)q" s -2a B {fluxcoeffb)q" b (3.65) 


11. Right-Front-Bottom Corner 

[a° p +®[2a E {coeffe) + a w +a N + 2a s {coeffs) + a T + 2a B {coeffb)^T P -®a w T w 
-© a N T N -®a T T T =(1 -®)[a w T* +a N T° +a T T T 0 ]+[a° p - {l-®)[2a E {coeffe) 

+ a w +a N + 2a s {coeffs) + a r + 2 a B {coeffb)]Tp 

+ [2a E {coeffe) +2a s {coeffs) + 2 a B (< zoeffb )JT cd -2a e { fluxcoeffe)q" e 

- 2a s {fluxcoeffs)q" s +2 a B {fluxcoeffb)q'\ (3.66) 


39 




12. Top Face 


\p p a s + 2a E (coefft) + a E jj/'p &a E T E 0a )( . T iV 

- 0 a N T N - 0a 5 r s - &a B T B = (1 - 0)[a £ r £ ° + a,T„° + a A ,7„ 0 + a s T s ° + a B T °] 

+ [a° P ~ (1 - ©)[^ £ + a w + a N +a s + 2 a T {coefft) + a £ ]]r £ 0 

+ [2a T (coefft)]r x - 2a T {fluxcoefft)cf\ -2 a T (fluxcoefft)q" S(mrce 

-2a r {fluxcoefft )q" boiling -2 a T (radcoeff)Tf (3.67) 


13. Front Face 

[a° p + 0[a £ + a w + a N + 2a s {coeffs) + a T +a B Jr, -Qa E T E -<da w T w 
- -®°tTt - ®a,T, = (1 -+ a w T£ + + a T 1? + a B T£] 

+ [ a l - (1 - ©)[ a £ + a w + a N + 2 a s (coeffs ) + a T +a B Jr; 

+ [ 2a s ( coeffs )jT co - 2 a s (fluxcoeffs)q ", (3.68) 


14. Back Face 

\flp { COe ff n ) a T Jr, — ® a E T E ~ ®a w T w 

-ea s T s -@a T T r -@a B T B =(l-0)[a £ r £ ° +a w T°+a s T s °+a r T r °+a B T°] 

+ [a° p - (1 - 0)[a £ + a w + 2 a N {coeffn) + a s + a T + a B Jr; 

+ [2a N {coeffn)]r x -2a N {fluxcoeffh)q\ (3.69) 


15. Left-Lateral-Face 

\a° p + 0[a £ + 2a w {coeffw) + a N + a s + a T + a B Jr, - Qa E T E - Qa N T N 
-&a s T s -0a r T r -&a B T B =(l-&)[aX+°X +°sT s ° +“ t T! +aj!] 

+ [a“ p -a- @)[a E + 2a u . (i coeffw ) + a N +a s +a r +a B Jr; 

+ \2a w {coeffwj\r x - 2a w (fluxcoeffw)q" w (3.70) 
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16. Right-Lateral-Face 


[a° p + ®[2a E (coeffe) + a w +a N +a s +a T + a B ^T P -®a w T w -®a N T N 
— ®a s T s — ®ci t T t — ®a B T B = (1 — + a N 7V + a s T s + a T T T + a B T B ] 

+ [< t ;-(l-®)[2a I ( coeffe ) + a w +a N +a s +a T +a B t? 

+ [2a E (coeffe)^ - 2a E (fluxcoeffe)q" e (3.71) 


17. Front-Top Edge 

[a° p + ©[a £ +a w + a N + 2 a s ( coeffs) + 2 a T (coeffi) + a B ^T P -®a E T E 

- ®a w T w — ®a N T N — ®a B T B = (1 — ©)[cz £ 7’ £ . + a w T w + o N T N + a B T B ] 

+ [a° p - (1 - ®)[a E +a w +a N + 2 a s ( coeffs ) + 2 a T (coeffi) + a B ^ff p 

+ [2a s ( coeffs ) + 2 a T (( coeffi )]T C0 - 2 a s (flvxcoeffs)q" s -2a T (fluxcoefft)q” t 

- 2a r (fluxcoefft)q" source -2a T (fluxcoefft)q\ oUing 

- 2a T (radcoeff)T p (3.72) 


18. Front-Bottom Edge 

[a° p + ®[a E +a w +a N + 2 a s (coeffs) + a T + 2a B (coeffbj§T p - ®a E T E 

— ®a w T w -®a N T N —®a T T T = (1 — +a w T w +a N T N + a T T T ] 

+ [a° p - (1 - ®)[a E +a w +a N + 2a s (coeffs) + a T + 2a B (coeffb)^T p 

+ [2« s (coeffs) + 2 a B (coeffb)]T x 

- 2 a s (fluxcoeffs)q" s -2a B (fluxcoeffb)q" b (3.73) 


19. Back-Top Edge 

[a° p + ®[a E +a w + 2a N (coeffn) + a s + 2a T (coeffi) + a B ^ff P -®a E T E 
-®a w T w - ®a s T s -®a B T B = (1 - ®)[a E T E +a w T° +a s T B + a B T B ] 

+ [a ° p (1 - ©)[tf £ + a w + 2a N (coeffn) + a s + 2a T (coeffi) + a B jr p ° 

+ \la N (coeffn) + 2a T (coeffi)^ -2a N (fluxcoeffh)q'\ -2 a T (fluxcoefft)q" t 
- 2a T (fluxcoefft )q" source ~2a r (fluxcoeffi)q" boiUng -2a r (radcoeff)T P (3.74) 
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20. Back-Bottom Edge 


[a° p + ©[fl £ + a w + 2a N {coeffri) + a s +a T + 2a B {coeffb)^T p - ®a E T E 
-®a w T w -®a s T s -Qa r T r = (1 -0)[a £ r £ ° + a w T* + a s T s ° + a T T T °] 
+ [ a ° P “ (1 “ ©)[«£ + a w + 2 a N {coeffri) + a s + a T + 2a B {coeffb)^T p 
+ [2a N {coeffri) + 2 a B ( coeffb)]T x 
- 2a N {flwccoeffn)q" n -2a B (fluxcoeffb)q" b 


21. Front-Right Edge 

[a° p + ©[2a £ ( coeffe ) + a w + a N + 2a s {coeffs ) + a T + a B jjr^ - ®a w T w 
-®a N T N -®a r T T -Ga B T B =(1 -®){a w T* +a N T° +a r T° + a B T B ] 
+ [a° p - (1 - ©)[2a £ {coeffe) + a w + a N + 2 a s {coeffs) + a T +a B ]]r £ 

+ [2a E {coeffe) + 2 a s {coeffs)]f x 
- 2a E {fluxcoeffe)q" e -2a s {fluxcoeffs)q" s 


22. Back-Left Edge 

[a° p + ®[a E + 2a w {coeffw) + 2a N {coeffri) + a s + a T + a B - ®a E T E 

- ®a s T s - ®a T T T - ®a B T B = (1 - ®)[a E T° + a s T s ° + a T T T ° + a B T° B ] 

+ [a° p - (1 - ©)[«£ + 2 a w {coeffw) + 2a N {coeffri) + a s +a T +a B ])t f 0 
+ [2a w {coeffw) + 2a N {coeffn)]T x 

- 2a w {fluxcoeffw)q'\ -2a N {fluxcoeffn)q" n 


23. Back-Right Edge 

[a° p +®[2a E {coeffe) + a w + 2 a N {coeffri) + a s + a T + a B f p -Ga w T w 
-®a s T s -®a r T T -®a B T B =(l-©)[a„T„° +a s Tf +a T Tf +a B T B ] 

+ l a ° P - (1 “ ®)\fa E {coeffe) + a w + 2a N {coeffri) + a s +a T +a B ~§T P 
+ [lo E {coeffe) + 2 a N {coeffn)]T x 
- 2 a E {fluxcoeffe)q" e -2a N {fluxcoeffn)q" n 


(3.75) 


(3.76) 


(3.77) 


(3.78) 
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24. Left-Lateral-Top Edge 

[a° p + ®\a E + 2 a w (coeffw) + a N +a s + 2 a T ( coefft) + a B lK - ®a E T E 

- ®a N T N - ®a s T s - Ga B T B = (1 - G)[a E T° + a N T° N + a s T s ° + affl ] 

+ [a° p - (1 - 0)[a £ + 2 a w (coeffw) + a N +a s +2a r (coefft ) + a B ]]r p 0 

+ [2 a w (coeffw) + 2 a T (coefft)}T x - 2a w (fluxcoeffw)q" w -2a T (fluxcoefft)q"\ 

- 2 a T (fluxcoefft)q” source -2a r (fluxcoejff)q" 

boiling 

- 2 a T (radcoeff)Tp (3.79) 


25. Left-Lateral-Bottom Edge 

[a° p + ®[a E + 2a w (coeffw) + a N + a s + a T + 2a B (coeffb)^T P -®a E T E 

- ®a N T N - ®a s T s - &a T T T = (1 - ®)[a E T° E + a N T* + a s T s ° + a T T T °] 

+ [a° p - (1 - ®)[a E + 2 a w (coeffw) + a N +a s +a T + 2a B (coeffb)\ (t p ° 

+ [2a w (coeffw) + 2 a B (coeffb)^ 

- 2a w (flvxcoeffw)q'\ -2 a B (fluxcoeffb)q'\ (3.80) 


26. Right-Lateral-Top Edge 

\a° p + ©[2 a E (coeffe) + a w +a N +a s + 2 a T (coefft) + a B lb - ®a w T w 

- ®a N T N - ®a s T s - ®a B T B = (1 - ®)[a w T* + a N T° N + a s T s ° + a B T° B ] 

+ [a° p - (1 - ©)[2 a E (coeffe) + a w + a N + a s + 2 a T (coefft) + a B Jr p ° 

+ [2a E (coeffe) + 2 a T (coefft)]T x - 2a E (fluxcoeffe)q" e -2a T (flvxcoefft)q", 

- 2 a T (fluxcoefft)q" source -2a T (fluxcoefft)q 

boiling 

-2a T (radcoeff)Tp (3-81) 
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27. Right-Lateral-Bottom Edge 


[a° p + ®[2a E (coeffe) + a w +a N +a s +a T + 2 a B (coeffb)\-®a w T w 
~ ®° N T N ~ ®a s T s - ®a T T T = (1 - ®)[a w T° + a N T° N + a s T s ° + a T T T ° ] 
+ [a° p - (1 - ©)[2 a E (coeffe) + a w + a N +a s +a T + 2 a B (coetfbjfr? 

+ [2a E (coeffe ) + 2 a B (coeffb)]T x 
- 2a E {fluxcoeffe)q" e -2a B (fluxcoeffb)q" b 


(3.82) 
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IV. RESULTS AND DISCUSSION 


Results are presented for the different cases discussed in greater detail below. As 
noted earlier, a variable sized and moving numerical mesh was used in such a way that 
the arc was always positioned at the center of the mesh where the spacing is the finest. 
The goal is to be able to resolve the large temperature gradient features around the arc 
and yet not incur a large overhead of computer resources, which would be required, if the 
grid was uniformly fine all over. This strategy however did require that a separate mesh 
generating routine be used which did demand some extra computational resources. The 
weld pool region was also modeled as a solid region but with a thermal conductivity 
higher than the surrounding unmelted region to simulate the effects of weld pool 
convection. The discontinuity in the thermal conductivity boundaries was handled using 
the standard technique of employing harmonic averaging at the boundary. Since the 
coefficients of the system of equations depend on the temperature, an iterative solution 
technique was used to achieve convergence in such a way that the maximum temperature 
difference between two consecutive iterations at any grid point was no more than 0.1 °C. 

The numerical solution method was used to examine different cases in freshwater 
for a 40-mm-thick 70 x 90 mm workpiece with a moving heat source in the positive y- 
direction. 

Case la: 

The validity of the numerical model was compared to Rosenthal’s three- 
dimensional solution for a moving heat source. At this point, convective, radiative and 
boiling surface thermal conditions were not considered. A constant thermal conductivity 
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value and a point heat source were used in the calculations. The input data used in 
computations for case la are shown in Table 1. The calculations were made up to 3.6 
seconds. The numerical results show excellent agreement with the analytical results of 
Rosenthal (Figure 4.1 and Figure 4.2). 


Workpiece : Length - 70 mm, Width = 90mm, Thickness = 40mm 
Water temperature: T = 27 °C (300.15 K) 

Power input into workpiece : Q = 2544 W 

Arc torch speed : v y = 4 mm/s 

Radius of heat input distribution : r 0 = 4.5 mm 

Thermal conductivity : k = 53 W/m K 

Density of steel: p = 7854 kg/m 3 

Specific heat of steel: C p = 509.3 J/kgK 


Table 1. The input data used in computations for case la 
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(b) 

Figure 4.2 The numerical model point heat source solution: 

(a) Profile; (b) Front-view. 

Case lb : 

In this case, the surface thermal boundary conditions of boiling heat transfer were 
now included. Other relevant data is given in Table 2. 

Workpiece : Length = 70 mm, Width = 90mm, Thickness = 40mm 
Water temperature: T = 27 °C (300.15 K) 

Saturation temperature: T=100°C (373.15 K) 

Water depth : 1 = Oft 

Total pressure: P = 101.325 kPa 

Power input into workpiece : Q - 2544 W 

Arc torch speed : v y = 4 mm/s 

Radius of heat input distribution : r 0 = 4.5 mm 

Thermal conductivity; k (W/m K) 

53 - 0.04 (T-300) 300 < T(K) < 1000 

25 + 6.25 x ia 3 (T-1000) 1000 < T(K) < 1800 

125 T(K) > 1800 

Emissivity : S - 0.82 _ 

Table 2. The input data used in computations for case lb 
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The top surface temperature values at 0.5 seconds can be seen in Figure 4.3 and 
Figure 4.4. At 0.5 seconds, the temperature distribution is almost symmetric about the 
position of the arc. The temperature distribution around the arc center is above the 
melting temperature of the workpiece (Tm = 1800 K). To see the depth of the weld pool 
penetration, the melting temperature contours were plotted for different surfaces (Figure 
4.5) and the weld pool depth was shown through the melting temperature points by using 
a curve-fit (Figure 4.6). Because of the brief reaction time, the arc heat input penetration 
to the workpiece is limited and a well-formed weld pool cannot be seen. [Ref. 16] 



Figure 4.3 Top surface t=0.5 sec 
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y-axis (m) 


Figure 4.4 Top surface (profile) t=0.5 sec 



Figure 4.5 The weld pool surface characteristics 
by the T=1800 K contours t=0.5 sec 
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Figure 4.6 The weld pool (front-view) t=0.5 sec 


At t = 2.0 sec., the results reveal that the temperature distribution has already 
reached steady state. Due to the moving heat source, the small change of slope of the 
temperature distribution curve can be observed clearly (Figure 4.7 and Figure 4.8). 
Because of the increased reaction time, the melting temperature contours are seen till the 
fourth surface from the top (Figure 4.9). The weld pool width increased to 8 mm and its 
depth to 2.1 mm (Figure 4.10). 
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Figure 4.7 Top surface t=2.0 sec 



Figure 4.8 Top surface (profile) t=2.0 sec 
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Figure 4.9 The weld pool surface characteristics 
by the T=1800 K contours t=2.0 sec 



-4 -3 -2 -1 0 12 3 4 

x-axis (mm) 


Figure 4.10 The weld pool (front-view) t=2.0 sec 
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At t= 4.5 sec., the change of slope of the temperature distribution behind the arc is 
more evident in Figure 4.11 and Figure 4.12. The melting temperature contours are still 
seen till the fourth surface from the top (Figure 4.13). But, the weld pool depth reached to 
3 mm with a width of 8 mm (Figure 4.14). These results show a good agreement with the 
data from Ref. 16. 



Figure 4.11 Top surface t=4.5 sec 
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Figure 4.13 The weld pool surface characteristics 
by the T=1800 K contours t=4.5 sec 
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Figure 4.14 The weld pool (front-view) t=4.5sec 


Case 2,3,4 : 

In the following cases we do not consider the temperature distribution. Instead, 

the cooling time that elapses between 800 °C and 500 °C is examined for freshwater at 
different temperatures and welding depths. The input data used in calculations are shown 
in Table 3, Table 4 and Table 5. 
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Workpiece : Length = 70 mm, Width = 90mm, Thickness = 40mm 
Water temperature: T=3°C (276.15 K) 

Saturation temperature : T = 114.88 °C (388.03 K) 

Water depth : 1 = 22 ft (6.706 m) 

Total Pressure : P = 168.75 kPa 

Power input into workpiece : Q = 3900 W 

Arc torch speed : v y = 2.75 mm/s 

Radius of heat input distribution : r 0 = 4.5 mm 

Thermal conductivity; k (W/m K) 

53 - 0.04 (T-300) 300 < T(K) < 1000 

25 + 6.25 x W 3 (T-1000) 1000 < T(K) < 1800 

125 T(K) > 1800 

Emissivity: £ = 0.82_ 


Table 3. The input data of calculations at T = 3 C and 1 = 22 ft 


Workpiece : Length = 70 mm, Width = 90mm, Thickness = 40mm 
Water temperature: T= 10 °C (283.15 K) 

Saturation temperature : T = 112.62 °C (385.77 K) 

Water depth : l = 18 ft (5.486 m) 

Total Pressure : P = 156.492 kPa 

Power input into workpiece : Q = 3900 W 

Arc torch speed : v y = 2.75 mm/s 

Radius of heat input distribution : r 0 = 4.5 mm 

Thermal conductivity; k (W/m K) 

53 - 0.04 (T-300) 300 < T(K) < 1000 

25 + 6.25 x W 3 (T-1000) 1000 < T(K) < 1800 

125 T(K) > 1800 

Emissivity: 8 = 0.82_ 


Table 4. The input data of calculations at T = 10 C and 1 = 18 ft. 
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Workpiece : Length = 70 mm, Width = 90mm, Thickness = 40mm 
Water temperature: T=31°C (304.15 K) 

Saturation temperature : T= 115.44 °C (388.59 K) 

Water depth : l = 24 ft (7.315 m) 

Total Pressure : P = 171.762 kPa 

Power input into workpiece : Q = 3900 W 

Arc torch speed : v y = 2.75 mm/s 

Radius of heat input distribution : r 0 = 4.5 mm 

Thermal conductivity ; k (W/m K) 

53 - 0.04 (T-300) 300 < T(K) < 1000 

25 + 6.25 x W 3 (T-1000) 1000 < T(K) < 1800 

125 T(K) > 1800 

Emissivity : 6 = 0.82_ 


Table 5. The input data of calculations at T = 31 C and 1 = 24 ft. 


The resulting graphs of the water temperature and welding depth effects on the 

peak temperature and the cooling time (between 800 °C and 500 °C ) for a point initially 
1.25 mm. behind the arc heat source are shown in Figure 4.15, Figure 4.16 and Figure 

4.17. The cooling times for 800 - 500 °C temperature range are 0.19 seconds(T wat er = 3 

°C, 1=22 ft.), 0.20 seconds (T wate r=10 °C, 1=18 ft.) and 0.33 seconds (T water =31 °C, 1=24 
ft.) respectively. These cooling times are significantly different from those in the previous 
studies in the literature. Based on the models in this case, it is seen that boiling 

phenomena must be considered for surface heat losses from the weld metal. For 800 °C 

to 500 °C range of the weld metal, film boiling exists and the surface is completely 
covered by a vapor blanket. A high amount of heat transfer from the weld metal to the 
water environment occurs by conduction across the vapor film. The other reason is that 

due to the thermal conductivity equations used in the model, for 800 °C to 500 °C 
ranges, thermal conductivity of steel is low. This causes a decrease in conduction heat 
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transfer from the arc to the studied point. But the high convection heat loss at this point 
causes a rapid cooling rate. Another possible reason is that because of the low welding 
speed and the high thermal conductivity in the weld pool, the workpiece, which has a 
large thermal capacity, absorbs most of the heat energy. After the arc has passed, 
temperature and thermal conductivity at that point drops. Due to the low thermal 
conductivity, the transfer of the absorbed heat energy from the inner part of the 
workpiece to the surface is not sufficient compared to the surface heat loss. This causes a 
rapid cooling rate on the surface. 

During welding, the high heat input to the workpiece increases the degree of grain 
growth. The grain coarsening effect causes a decrease in the grain boundary area. 
Because of the large grain size and the very high cooling rates, the resulting phase is hard 
with Vi > 450, but with a brittle martensitic structure. 



Figure 4.15 Cooling time and peak temperature for a point 
initially 1.25 mm behind the arc heat source. 
(T water = 3 C, depth = 22 ft.) 
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Figure 4.16 Cooling time and peak temperature for a point 
initially 1.25 mm behind the arc heat source. 
(Twater ~ 10 C ( depth = 18 ft.) 



Figure 4.17 Cooling time and peak temperature for a point 
initially 1.25 mm behind the arc heat source. 
(T W ater = 31 C, depth = 24 ft.) 
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V. CONCLUSIONS AND RECOMMENDATIONS 


Compared to experimental methods, numerical solution techniques are very fast, 
effective and economical alternatives. Therefore, the aim of this study was to develop a 
general numerical model for transient, three-dimensional conduction heat transfer 
phenomena in underwater welding process on a thick rectangular plate. Computations 
were presented for different cases. Comparisons of current predictions with results in the 
literature showed good agreement and validated the model. If the material, environment 
and the arc source properties are known, this program can be applied to the different 
types of metals under the wet welding process. 

The numerical model gave important temperature-time data in the critical HAZ, 
which in turn determines material structure. The model also helped to make prediction of 
size of weld pool and its evolution with time. 

The principal recommendations for the future studies may be summarized as 
follows. 

1. The weld pool region itself was modeled, as a solid region. But, owing to 
the presence of the temperature values higher than the melting temperature 
in the weld pool, the model must be considered as a liquid weld pool with 
the surrounding unmelted solid region. So, the numerical model needs to 
be modified to account for weld pool convection with the buoyancy force, 
the electromagnetic force and the surface tension gradient at the weld pool 
surface. 

2. The chemical effects due to interaction of the arc with the surrounding 
water environment may need to be included. 

3. In the model, boiling heat transfer phenomena was accounted for a solid 
surface. For a molten pool region, the role of boiling on a liquid substrate 
is unknown and needs to be solved. 
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APPENDIX A. PROGRAM STRUCTURE 


In the main program <aathl0.m>, non-uniform grid spacing was constructed by 
using the exponential function y = he~ bx where b is the coefficient which affects the grid 
spacing distribution and h is the sum of the grid spacings when x -»oo. The application 
of grid spacing to the workpiece was shown in Figure Al. 



(b) 

Figure Al Applying the grid spacing to the workpiece by using 
the exponential function 
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The distance in front of the moving heat source was called as “front” and the 
distance behind the heat source was called as “back”. The width and the thickness of the 
plate were defined with the same terms. The variable “number” was used to define the 
number of grid points in the “back” region. Variables “b” and “number” can be changed 
to control the resolution of the grid spacing. 

To simulate the heat input from the arc to the top surface, it was assumed that the 
arc source heat flux distribution had a Gaussian distribution on the top surface. The 
applied arc heat source was defined as a matrix, which its size was equal to the size of the 
top surface. The value of heat flux at each grid point was calculated by using the x- 
coordinate and the y-coordinate of that point. The resulting arc heat source matrix was 
applied to the top surface. 

The discretised equations was represented by the matrix equation [A] X =b 
where, [A] is the coefficient matrix, X is the column matrix of the unknown temperature 
values of the grid points and b is the column matrix of the constants. In the coefficient 
matrix, the coefficients of the studied grid and the neighboring grids were written to the 
same row. An example for 5x5x2 workpiece was shown in Figure A2. In this example, 
the coefficient matrix may be written as 

cl c2 0 0 0 c6 0... c26... 

cl c2 c3 0 0 0 c7 0... c27... 

0 c2 c3 c4 0 0 0 0... c28... 

0 0... c25... c45... c49 c50 
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Figure A2 The example workpiece 


Even the row number, the column number and the surface number were increased, 
the m ax im um number of coefficients for each row does not exceed 7. This is for an 
interior grid point of a workpiece with a surface number 3 or more. But, the increasing 
size of the coefficient matrix increases the computing time and the required computing 
memory and storage capacity. To prevent this problem, tridiagonal part of the general 
matrix (Figure A3) which contains T w , T P and T E coefficients was formed as a new ax3 

matrix (a = row x column x surface). The coefficient values and T n , T s , T t and T B 

were also formed separately and called as <north>, <south>, <top> and <bottom> 
matrices respectively. The size of each matrix is a x 1 (a = row x column x surface). The 
constant values in the right side of the equation were called as <right> matrix. 
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Figure A3 The coefficient matrix 


Beginning from the first surface, the function program <aasolvelO> solves the 
matrix produced by the main program by using row-by-row iterative sweeping method. 
The representation of row-by-row method was shown in Figure A4 [Ref. 21,22]. The 
statements <indexl> and <index2> defines the first and the last element of each row. 



• Points at which values are calculated 

O Points at which values are considered 
to be temporarily known 

■ Known boundary values 
—► Sweeping direction 

Figure A4 Row-by-row method representation 






The matrix <temp_old> represents the old temperature values of the grid points. 
At first, it was initialized and then it took the values of the previous iteration. The matrix 
<temp> represents the temperature values of the grid points, which were found from the 
present iteration. The sizes of <temp_old> and <temp> are (a x b) x c, where a, b and c 
are the total row number, column number and surface number respectively. The statement 
<delta_iteration> represents the absolute value of the difference between <temp_old> and 
<temp>. If <delta_iteration> is less 0.1, the resulting temperature value is satisfactory. 
The statement <count_iteration> counts the iteration number. If the program can not 
converge after 10 iterations, <delta_iteration> increases 0.1 and continues to increase 0.1 
at every following 5 iterations until the program find a converged temperature value. The 
radiation heat flux matrix <radiationl> and the boiling heat flux matrix <q_boiling_free> 
are only used to calculate the heat flux from the top surface.. The sizes of <radiationl> 
and <q_boiling_free> are (a x b) x 1, where a, b are the total row number and column 
number respectively. 

In the <aasolvel0> function program, <north>, <south>, <top> and <bottom> 
matrices were passed to the right side of the equation. Their values were calculated by 
using the temperature values from the previous iterations and the resulting values were 
added to the values of <right>. The diagonal matrix was solved by using Gauss 
elimination method and the values of T w , T P and T E was found for the present iteration. 

In the <old_templ0> function program, it was assumed that the plane was moving 
with the arc source together (Figure A5). Here, due to the non-uniform grid spacing, it 
was necessary to determine the new position and the temperature value of each grid point 
by using second degree polynomial enterpolation (Figure A6). The values of T a , T b and 
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T c were found by using the old temperatures from the previous time step. These 
temperatures was used to determine the new grid temperature T x with the second degree 
polynomial enterpolation. 
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Figure A5 Moving plane 



(a) 


(b) 


Figure A6 Grid positioning by using second degree polynomial 
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The function program <save_alllO> saves all 3-D temperature data at each time 
step by overwriting onto. The other function program <save_thelO> saves top surface 
temperatures at each time step into a different file name. 
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APPENDIX B. PROGRAM CODES 


THE MAIN PROGRAM (aathlO.m) 


% Need to define the function below to be able to use mcc 
% The symbol "al" itself has no significance - aathlO is the 
% name of the main program, 
function [al]=aathlO () 

clear % clear all variables - not present in "aathlOgoon .m" 

close all % close all open figure windows 

% save_me is a counter that is incremented after each time step 
% and is used to decide how often to save the calc temperature 
% data. 
save_me=0; 

% save_counter is a counter that is also incremented with each 
% time step and is appended to the file name into which the temp 
% data is stored. 

save_counter=0; % not present in "aathlOgoon.m" 

% Tinf is the ambient temp in deg C 
Tinf=27; 

% TO is the initial temp in deg C 
T0=2 7; 

% sil is the stefan-boltzmann const (SI units) 
sil=5.67E-8; 

% epsilon is the emissivity of the surface 
epsilon=0; 

% x-vel component of torch speed (m/s) 
vel_x=0; 

% y-vel component of torch speed (m/s) 
vel_y=0.004; 

% deltat_t is the time step (in secs) 
delta_t=0.001; 

% Distance source moves in the x-direction in each time step (in m) 
deltax=vel_x*delta__t ; 


% Distance source moves in the y-direction in each time step (in m) 
deltay=vel__y*delta_t ; 

% q_up is the heat flux distribution on the upper surface (W/m^2) 
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% q_up is the heat flux distribution on the upper surface (W/m A 2) 
% imposed by the torch/arc (implemented in matrix form) 
q_up=0; 


%% q is a generic variable representing heat flux that may be 
%% present at any of the other surfaces and can be included 
%% in the nodal equations 

% q_ is constant heat flux through the surfaces 

q_west=0; 

q_east=0; 

q_top=0; 

q_bottom=0; 

q_north=0; 

q_south=0; 

% c is the cp [J/kg-K] 
c-509.3; 

% ro is the density [kg/m A 3] 
ro=7854; 

% Convection heat transfer coefficients from the faces [W/m A 2-K] 

he=0; 

hw=0; 

hn=0; 

hs=0 ; 

ht=0; 

hb-0; 

% The distances between the grid point points 

%% back is the portion of the moving grid behind the point 

%% of location of the source [in m] 

back=.05; 

% b is the coefft in the exponential gridding function (e A (-b*x)) 
b=.13; 

% number is the no. of grid points in the "back” region 
number=20; 


%% Variables "b" and "number" above can be varied to control the 
%% features of the variable grid, such as change in grid spacing, 

%% etc. 

%% height is a vector that holds the coords of the grid points in %% 
the back-region. Note that the spacing between these grid 
%% points is the diff. between consecutive height entries and is 
%% stored in "int_back" 
height(1)=0; 
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% Note that int__back(l) > int_Joack(number) , i.e. in reverse order 
for kk=l:number 

height(kk+1)=back*(1-exp(-b*kk) ) ; 
int_back(kk)=height(kk+1)-height(kk); 

end 


%% multiply is a scaling factor(close to 1) that scales height and 
%% int_back in such a way that it height(number) = back. 
multiply=back/(sum(int_back)-intjoack(l) 12 ) 
int_back=multiply*int_back; 

%% Similar to back, below are the variables front (front region of %% 

%% moving grid ahead of source), width (width of plate) and 

%% thickness. All distances in m. 

front=0.02; 

width=0.09; 

thickness=0.04; 

%% int_front is the grid spacing in the front region. Initially 

%% set to 0. Then transfer from int-back one by one until sum of 

%% int-front entries exceeds front. 

int_front=0; 

counter=0; 

while (sum(int_front)-int_front(length(int_front))/2)<=front 


% int_front below continually changes size as the spacings are 
% added. 

int_front=[int_front,int_back(length(int_back)-counter)]; 
counter=counter+l; 

end 

% Set first member which was initially 0 to null. 

% Meaningful spacing starts only from the 2nd member of int_front 
int front(!)=[]; 


%% The grid distance matrix in the y direction is y_intl starting 
%% from the top of front down through the source all the way to 
%% the bottom of back. 

% fliplr is to ge the desired ordering from front to back. 
y_intl=[fliplr(int_front),fliplr(int_back)]; 

% int_side holds grid spacing values over the width. 

% Same logic as for int_front 

int_side=0; 

counter=0; 

while (sum(int_side) -int__side (length (int__side) ) /2) <= (width/2) 
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int__side= [int_side, int_back (length (int_back) - counter) ] ; 
counter=counter+l; 
end 

int_side(!)-[]; 


% The grid distance vector in the x dirextion is x_intl (y_intl 
% above. 

x_intl=[fliplr(int_side),int_side]; 

% The grid distance vector in the z dirextion is z_intl (y_intl 

% above) 

z_intl=0; 

counter=0; 

while (sum (z_intl) -z__intl (length (z_intl) ) /2) <=thickness 

z_intl=[z_intl,int_back(length(int_back)-counter)]; 
counter=counter+l; 


end 

z_intl(!)=(]; 


figure(1) 

subplot (3,1,1) 
plot(x_intl) 
grid on 


subplot (3,1,2) 
plot(y_intl) 
grid on 

subplot(3,1,3) 
plot(z_intl) 
grid on 

% row is the number of y nodes (along the length). 
row = length (y__intl)-1 ; 

% col is the no. of x-nodes (across the plate width), 
col = length(x_intl)-1; 

% surface is the no. of z-nodes (across the thickness). 
surface = length(z_intl)-l; 
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The value of teta determines the numerical scheme as follows: 


% for Fully Implicit Method; teta=l 
% for Crank-Nicholson Method; teta=0.5 
% for Explicit Method; teta=0 

teta=l; 


For a a*b*c three dimensional matrix 


% a:num rows, b:num cols, c:num 


% 

o. 

Size 

of 

new- 

—> 

(a*b*c)*3 

*0 

% 

Size 

of 

top- 

—> 

(a*b*c)*1 

% 

% 

size 

of 

bottom 

—> 

(a*b*c)*1 

% 

% 

size 

of 

north - 

—> 

(a*b*c)*1 

% 

% 

size 

of 

south - 

— > 

(a*b*c)*1 

% 

% 

a 

size 

of 

right - 

— > 

(a*b*c)*1 

o 

% 

size 

of 

temp — 

—> 

(a*b)*c 


surfaces. 

TP, TE, TW coefficients for the 
whole grids. 

TT coefficients for the whole 
grids. 

TB coefficients for the 
whole grids. 

TN coefficients for the whole 
grids. 

TS coefficients for the whole 
grids. 

The values of the righthandside 
of the matrix. 

The temperatures of the grids. 


% new is the coefft matrix for all the grid points in the 3-D 
% volume. 


% 1st col holds the coeffts of the West (W) nodes. 

% 2nd col holds the coeffts of the (P) nodes. 

% 3rd col holds the coeffts of the East (E) nodes. 

new=zeros(row*col*surface, 3); 

% right is the const (or b-vector) on the RHS 
right=zeros(row*col*surface,1); 

%% south, north, top, bottom hold the coeffts of the corresponding %% 

grid points. 

south=right; 

north=right; 

top-right; 

bottom=right; 

% Apply the initial condition to the plate. 

% "temp" is the temperature matrix. 

% temp is a matrix that holds temperature values for each surface. 

% Note dimensions of temp carefully. 

% This initialization is not present "aathlOgoon.m". 
temp=T0*ones(row*surface,col); 
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% Below are parameters that specify the heat/arc source. 

% rO is the effective radius [m] of the gaussian heat source. 
r0=0.0045; 

% ddd is a coefft that determines the gaussian spread. 

% The larger the ddd value, the more pointed the source. 
ddd=3; 

% qrate_total is the total energy transfer rate from the source 
% [Watts]. 
qrate_total = 2544; 


% qO is the peak heat flux of the gaussian [W/m A 2] 
qO=qrate__total*ddd/pi/rO A 2; 

% Calculate the q_up (heat flux distribution on the top face) 

%% xcenter and ydenter are the grid element numbers of the grid 
%% center at which the source/arc is located. 
xcenter=length(int_side); 
ycenter=length(int_front); 


for nn=l:row 
for mm=l:col 

% x(mm) is array with the x-coords of the mesh, 
if mm < xcenter 

x(mm)=-sum(x_intl((mm+1):xcenter)); 
elseif mm==xcenter 
x (mm) =0; 
else 

x(mm)=sum(x_intl({xcenter*1):mm)); 

end 

% y(mm) is array with the y-coords of the mesh, 
if nn < ycenter 

y(nn)=sum(y_intl((nn+1):ycenter)); 
elseif nn—ycenter 
y(nn)=0; 
else 

y(nn)=-sum(y_intl((ycenter+1):nn)); 

end 

% initialize q_up using x() and y() coord info. 

q_up(nn,mm)=-q0*exp(-ddd*(x(mm) A 2+y(nn) A 2)/rO A 2); 

end 

end 

figure 

% a mesh command to view the generated grid, 
mesh (x, y, q_up) 
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for time=delta_t:delta_t:200 

% 200 is the final time in secs upto which the program should 
% go. . . 

save_me=save_me+l; 

disp('the interpolation time is 1 ) 
tic % beginning of tic-toc loop 

[temp]=old_templ0(row,col,surface,temp,x,y,deltax,deltay); 

toe % end of tic-toc loop. 


% figure 

% contour(x,y,temp(1:row,1:col) ) 
disp('The matrix formation time is : ? ) 
tic % beginning of tic-toc loop. 

%% Start to form the matrix by defining every point. 

% Beginning of creation of coefft matrix 

% loops step through each point on all the surfaces, cols and rows 

for kk=l:surface 
for jj=l:row 
for ii=l:col 


% The distances between the reference point and the adjoining 
% nodal points. 

% Distance between P&E and P&W 
xpe=x_intl(ii+1); 
xwp=x_intl(ii) ; 

% Distance between P&N and P&S 
ypn=y_intl{jj); 
ysp=y_intl(jj+1); 

% Distance between P&T and P&B * 
zpt=z_intl(kk); 
zbp=z_intl(kk+1) ; 

%% The lateral surface areas of the faces of the control volume 
%% around the reference point.P. 

% East-face area. 

Ae=(ypn/2+ysp/2 )* (zpt/2+zbp/2); 

Aw=Ae; 
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% North-face area. 

An=(xpe/2+xwp/2)*(zpt/2+zbp/2); 

As=An; 

% Top-face area. 

At=(xpe/2+xwp/2)*(ypn/2+ysp/2) ; 

Ab-At; 

% The size of the control volume. 

deltaV=(xpe/2+xwp/2)*(ypn/2+ysp/2)*(zpt/2+zbp/2); 


% Routine to set thermal conductivities below: 

% Initialize TEMP ()s to 0 for logical ifs below: 

TEMP(1)=0;TEMP(2)=0;TEMP(3)=0;TEMP(4)—0;TEMP(5)=0;TEMP(6)=0; 


% TP0 is the temperature of the current node. 

TP0=temp{(kk-1)*row+jj, ii) ; 

% 

% 

% T_west=TEMP(1) T_east=TEMP(2) T_north=TEMP(3) 

% T_south=TEMP(4) T_top=TEMP (5) T_bottom=TEMP ( 6) 

% 

% 

if ii==l % first col, no west neighbor, so... 

TEMP(1)=TP0; % set T_west to TP0 

elseif ii==col % last col, no east neighbor, so. 

TEMP (2 ) =TP0; % set Tjast to TP0 

end 

if jj—1 

TEMP(3)=TP0; 
elseif jj==row 
TEMP(4)=TP0; 

end 

if kk==l 

TEMP(5)=TP0; 
elseif kk==surface 
TEMP(6)=TP0/ 

end 

%% If current node P has a valid neighbor, then use temperature 
%% of that node.... 

if (TEMP(1)==0) 

TEMP(1)=temp((kk-1)*row+jj , ii-1) ; 

end 
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if (TEMP(2)==0) 

TEMP(2)=temp((kk-1)*row+jj , ii+1) ; 

end 

if (TEMP(3)==0) 

TEMP(3)=temp((kk-1)*row+jj-1,ii); 
end 

if (TEMP(4)==0) 

TEMP(4)=temp((kk-1)*row+jj + 1, ii); 
end 

if (TEMP(5)==0) 

TEMP(5)=temp((kk-2)*row+jj , ii); 
end 

if (TEMP(6)==0) 

TEMP(6)=temp(kk*row+jj ,ii); 
end 


% Do loop to calculate k-values of all grid points. 

for kjl=l:6 

% Now set the thermal conductivity values. 
kelvin=TP0+273.15; 

kl=53; % nominal kl-value in W/m-K set here, but... 

% Actual kl-values are calc below, 
if (kelvin>=300 & kelvin<1000) 
kl=53-0.04*(kelvin-300); 
elseif (kelvin>=1000 & kelvin<1800) 
kl=25+6.25E-3*(kelvin-1000); 
elseif kelvin>=1800 
kl-125; 

end 

% Find the thermal conductivity at the adjacent point 
kelvin=TEMP(kj1)+273.15; 

k2=53; % nominal k2-value in W/m-K set here, but... 

% Actual k2-values of neighbors calc here, 
if (kelvin>=300 & kelvin<1000) 
k2=53-0.04*(kelvin-300) ; 
elseif (kelvin>=1000 & kelvin<1800) 
k2=25+6.25E-3*(kelvin-1000) ; 
elseif kelvin>=1800 
k2=125; 

end 
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% Harmonic mean thermal conductivity found next 
NEW__K(kjl) = 2*kl*k2/ (kl + k2) ; 

end 

% Now redefine k f s of the neighboring faces with values from NEW K 

kw=NEW_K(1); 

ke=NEW_K(2); 

kn=NEW_K(3); 

ks=NEW_K(4) ; 

kt==NEW_K (5) ; 

kb=NEW__K(6) ; 

% Coefficients used in developing the discretised equations. 

ae=ke*Ae/xpe; 

aw=kw*Aw/xwp; 

an=kn*An/ypn; 

as=ks*As/ysp; 

at=kt*At/zpt; 

ab=kb*Ab/zbp; 

apO=ro*c*deltaV/delta_t; 

ap=teta* (ae+aw+an+as+at+ab) + apO; 

% Coefficients used for the wall-medium interface part of the 

% discretised equations. 

coeffe=he*xpe/(2*ke+he*xpe); 

coeffw=hw*xwp/(2*kw+hw*xwp); 

coeffn=hn*ypn/(2*kn+hn*ypn); 

coeffs=hs*ysp/(2*ks+hs*ysp); 

coefft=ht*zpt/(2*kt+ht*zpt); 

coeffb=hb*zbp/(2*kb+hb*zbp); 

% Coefficients used for the heat flux part of the discretised 
% equations. 

fluxcoeffe=xpe/(2*ke+he*xpe); 
fluxcoeffw=xwp/(2*kw+hw*xwp); 
fluxcoeffn=ypn/(2*kn+hn*ypn); 
fluxcoeffs=ysp/(2*ks+hs*ysp); 
fluxcoefft=zpt/(2*kt+ht*zpt); 
fluxcoeffb=zbp/(2*kb+hb*zbp); 

% Radiation coefft for surface b.c. - see thesis... 
radcoeff=sil*epsilon*zpt/(2*kt+ht*zpt); 


% row_num is the index of the current point in the matrices being % 
used. row_num varies from 1:surface*col*row 
row__num= (kk-1) *row*col+ (j j -1) *col+ii; 

% For left-back-top corner 


if ((kk==l) & (jj=l) & (ii==l) ) 
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% Temperature values for current time step for each relevant node. 
TP0=temp((kk-1)*row+jj , ii) ; 

TE0=temp((kk-1)*row+jj , ii+1) ; 

TS0=temp ( (kk-1)*row+j j+1, ii) ; 

TB0=temp(kk*row+jj , ii) ; 

% Coeffts of grid points; note south and bottom later are moved to 
% the rhs. ■ 

% "new" is the A-matrix in the solution. 

% Note teta is the parameter which determines the numerical method 
% teta=l:fully implicit; teta=0:explicit; teta=0.5:Crank-Nicholson 
new(row_num,2)=apO+teta*(ae+2*aw*coeffw+2*an*coeffn+as+2*at*coefft+ab) 
new(row_num,3)=-teta*ae; 
south(row_num,1)=teta*as; 
bottom(row_num,1)=teta*ab; 

% "right" is the constant b-matrix on the rhs... 
right(row_num,1)=(1-teta)*(ae*TE0+as*TS0+ab*TB0)+ ... 

(apO-(1-teta)*(ae+2*aw*coeffw+2*an*coeffn+as+2*at*coefft+ab))*TP0+ .. 
(2*aw*coeffw+2*an*coeffn+2*at*coefft)*Tinf - 
2*at*fluxcoefft*q_up(j j , ii) ... 

-2*aw*fluxcoeffw*q_west-2*an*fluxcoeffn*q_north-2*at*fluxcoefft*q_top 

%% radiation & q_boiling_coeff belong to the rhs but are not included 
%% until later since they depend on the current temperature 


% (a form of quasi-linearization of nonlinear temp terms) 
radiation(row_num,1)=-2*at*radcoeff; 
q_boiling_coeff(row_num,1)=-2*at*fluxcoefft; 

% for right-back-top corner 

elseif ((kk—1) & (jj==l) & (ii==col) ) 

TP0=temp((kk-1)*row+jj, ii) ; 

TW0=temp((kk-1)*row+jj , ii-1) ; 

TS0=temp((kk-1)*row+jj+1,ii) ; 

TB0=temp(kk*row+jj,ii); 

new (row_num, 2)=apO+teta* (2*ae*coeffe+aw+2*an*coeffn+as+2*at*coefft+ab) 
new (row__num, 1) =-teta*aw; 
south(row_num,1)=teta*as; 
bottom(row_num,1)=teta*ab; 

right(row_num,l)=(l-teta)*(aw*TW0+as*TS0+ab*TB0)+ ... 

(apO-(1-teta)* (2*ae*coeffe+aw+2*an*coeffn+as+2*at*coefft+ab))*TP0+ 

(2*ae*coeffe+2*an*coeffn+2*at*coefft)*Tinf- 
2*at*fluxcoefft*q_up(jj , ii). . . 

-2*ae*f luxcoef fe*q__eas t-2* an* f luxcoef fn*q_north- 


2*at*fluxcoefft*q_top; 

radiation (row_num, 1) =-2*at*radcoeff; 

. 81 




q_boiling_coeff (row__num, 1) =-2*at*fluxcoef ft ; 

% for left-front-top corner 

elseif ((kk~=l) & (jj==row) & (ii—=1)) 

TP0=temp((kk-1)*row+jj,ii); 

TE0=temp((kk-1)*row+jj , ii+1); 

TN0=temp((kk-1)*row+jj-1, ii); 

TB0=temp(kk*row+jj,ii); 

new(row^num,2)=apO+teta*(ae+2*aw*coeffw+an+2*as*coeffs+2*at*coefft+ab); 
new(row_num,3)=-teta*ae; 
north(row_num,1)=teta*an; 
bottom(row_num,1)=teta*ab; 

right(row_num,1)=(1-teta)*(ae*TE0+an*TN0+ab*TB0)+ ... 

(apO-(1-teta)*(ae+2*aw*coeffw+an+2*as*coeffs+2*at*coefft+ab))*TP0 + ... 
(2*aw*coeffw+2*as*coeffs+2*at*coefft)*Tinf-2*at*fluxcoefft*q_up(jj,ii) 

-2*aw*fluxcoef f w*q__west-2*as* f luxcoef fs*q_south-2*at*f luxcoef ft *q top; 
radiation(row_num,1)=-2*at*radcoeff; ™ 

q_boiling_coeff(row__num,1)=-2*at*fluxcoefft; 


% for right-front-top corner 

elseif ((kk==l) & (jj==row) & (ii==col)) 

TP0=temp((kk-1)*row+j j , ii) ; 

TW0=temp((kk-1)*row+j j , ii-1); 

TN0=temp{(kk-1)*row+jj-1, ii) ; 

TB0=temp(kk*row+j j , ii) ; 


new(row_num,2)=apO+teta*(2*ae*coeffe+aw+an+2*as*coeffs+2*at*coefft+ab); 
new(row_num, 1) =-teta*aw; 
north(row_num,1)=teta*an; 
bottom (row__nura, 1) =teta*ab; 

right(row^num,1)=(1-teta)*(aw*TW0+an*TN0+ab*TB0)+ ... 

(apO-(1-teta)*(2*ae*coeffe+aw+an+2*as*coeffs+2*at*coefft+ab))*TP0 + 

(2*ae*coeffe +2*as*coeffs +2*at*coefft)*Tinf- 
2*at*fluxcoef ft *q__up {j j , ii) ... . 

-2*ae*fluxcoeff e*q_east-2*as*fluxcoeffs*q_south- 
2*at *fluxcoef f t*q_top; 
radiation(row_num,1)=-2*at*radcoeff; 
q__boiling_coef f (row^num, 1) =-2 *at *f luxcoef ft; 
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% for left-back-bottom corner 


elseif ((kk==surface) & (jj==l) & (ii==l)) 

TP0=temp((kk-1)*row+jj , ii) ; 

TE0=temp((kk-1)*row+jj,ii+1); 

TS0=temp((kk-1)*row+jj+1,ii) ; 

TT0=temp((kk-2)*row+jj, ii); 

new(row_num,2)=apO+teta*(ae+2*aw*coeffw+2*an*coeffn+as+at+2*ab*coeffb) 
new(row_num,3)=-teta*ae; 
south(row_num,1)=teta*as; 
top(row_num,1)=teta*at; 

right(row_num,1)=(1-teta)*(ae*TE0+as*TS0+at*TT0)+ ... 

(apO-(1-teta)*(ae+2*aw*coeffw+2*an*coeffn+as+at+2*ab*coeffb))*TP0+ .. 
(2*aw*coeffw+2*an*coeffn+2*ab*coeffb)*Tinf ... 
-2*aw*fluxcoeffw*q_west-2*an*fluxcoeffn*q_north- 
2*ab*fluxcoeffb*q_bottom; 


% for right-back-bottom .corner 

elseif ((kk==surface) & (jj==l) & (ii==col)) 

TPO^temp((kk-1)*row+jj , ii) ; 

TW0=temp((kk-1)*row+jj,ii-1); 

TS0=temp((kk-1)*row+jj + 1,ii) ; 

TTO-temp((kk-2)*row+jj , ii) ; 

new(row_num,2)=apO+teta*(2*ae*coeffe+aw+2*an*coeffn+as+at+2*ab*coeffb) 
new(row__num, 1) =-teta*aw; 
south(row_num r 1)=teta*as; 
top(row_num,1)=teta*at; 

right(row_num,1)-(1-teta)*(aw*TW0+as*TS0+at*TT0)+ ... 

(apO-(1-teta)*(2*ae*coeffe+aw+2*an*coeffn+as+at+2*ab*coeffb))*TP0+ 

(2*ae*coeffe+2*an*coeffn+2*ab*coeffb)*Tinf ... 

-2*ae*fluxcoeffe*q_east-2*an*fluxcoeffn*q_north- 
2*ab*fluxcoeffb*q_bottom; 


% for left-front-bottom corner 

elseif ((kk==surface) & (jj==row) & (ii==l)) 

TPO-temp((kk-1)*row+jj , ii) ; 

TE0=temp((kk-1)*row+j j,ii + 1); 

TN0=temp((kk-1)*row+jj-1,ii); 

TT0=temp((kk-2)*row+jj , ii) ; 

new(row_num,2)=apO+teta*(ae+2*aw*coeffw+an+2*as*coeffs+at+2*ab*coeffb) 
new(row num,3)=-teta*ae; 
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north(row_num,1)=teta*an; 
top(row_num,1)=teta*at; 

right(row_num,1)=(1-teta)*(ae*TE0+an*TN0+at*TT0)+ . . . 

(apO-(1-teta)*(ae+2*aw*coeffw+an+2*as*coeffs+at+2*ab*coeffb))*TP0+ .. 
(2*aw*coeffw+2*as*coeffs+2*ab*coeffb)*Tinf ... 

-2*aw*fluxcoeffw*q_west-2*as*fluxcoeffs*q_south- 
2*ab*fluxcoeffb*q_bottom; 


% for right-front-bottom corner 

elseif ((kk==surface) & (jj==row) & (ii==col)) 

TP0=temp((kk-1)*row+j j ,ii); 

TW0=temp((kk-1)*row+jj,ii—1) ; 

TN0=temp((kk-1)*row+jj-1, ii); 

TT0=temp((kk-2)*row+jj,ii); 

new(row_num,2)-apO+teta*(2*ae*coeffe+aw+an+2*as*coeffs+at+2*ab*coeffb) 
new(row__num, 1) =-teta*aw; 
north (row_num, 1) =teta*an; 
top(row^num,1)=teta*at; 

right(row_num,1)=(1-teta)*(aw*TW0+an*TN0+at*TT0)+ ... 

(apO-(1-teta)*(2*ae*coeffe+aw+an+2*as*coeffs+at+2*ab*coeffb))*TP0+ 

(2*ae*coeffe +2*as*coeffs +2*ab*coeffb)*Tinf ... 

-2*ae*fluxcoeffe*q_east-2*as*fluxcoeffs*q_south- 
2 *ab*f luxcoef fb*q_bottom; 


% for the top surface 

elseif ((kk==l) & (j j >1) & (jj<row) & (ii>l) & (ii<col)) 

TP0=temp((kk-1)*row+j j , i i) ; 

TE0=temp((kk-1)*row+jj,ii+1); 

TW0=temp((kk-1)*row+jj,ii-1); 

TN0=temp((kk-1)*row+jj-1,ii); 

TS0=temp((kk-1)*row+jj+1,ii); 

TB0=temp(kk*row+jj, ii) ; 

new(row_num,2)=apO+teta*(ae+aw+an+as+2*at*coefft+ab) ; 

new(row_num,3)=-teta*ae; 

new(row_num,1)=-teta*aw; 

north(row_num,1)=teta*an; 

south (row_num, l)=teta*as; 

bottom(row_num,1)=teta*ab; 

right(row_num,1)-(1-teta)*(ae*TE0+aw*TW0+an*TN0+as*TS0+ab*TB0)+ ... 
(apO-(1-teta)*(ae+aw+an+as+2*at*coefft+ab))*TP0 + 2*at*coefft*Tinf - 

2*at*fluxcoefft*q_up(jj,ii)-2*at*fluxcoefft*q_top; 
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radiation(row_num,1)=-2*at*radcoeff; 
q_boiling_coeff(rowjium,1)=-2*at*fluxcoefft ; 


% for the bottom surface 

elseif ((kk==surface) & (jj>1) & (jj<row) & (ii>l) & (ii<col)) 

TP0=temp ( (kk-l)*row+jj , ii) ; 

TE0=temp((kk-1)*row+jj,ii+1); 

TW0=temp((kk-1)*row+jj,ii-1) ; 

TN0=temp((kk-1)*row+jj-1, ii); 

TS0=temp((kk-1)*row+jj+1,ii) ; 

TTO^temp((kk-2)*row+jj, ii); 

new (row_num, 2)=apO+teta*(ae+aw+an+as+at+2*ab*coeffb) ; 

new(row_num,3)=-teta*ae; 

new (row__num, 1)=~teta*aw; 

north(row_num,1)=teta*an; 

south (row__num, 1) =teta*as; 

top (row__num, 1) =teta*at; 

right(row^num,1)=(1-teta)*(ae*TEO+aw*TWO+an*TNO+as*TSO+at*TTO)+ ... 
(apO-(1-teta)*(ae+aw+an+as+at+2*ab*coeffb))*TP0 + 2*ab*coeffb*Tinf 

2*ab*fluxcoeffb*q_bottom; 


% for the front surface 

elseif ((kk>l) & (kk<surface) & (jj==row) & (ii>l) & (ii<col)) 

TP0=temp((kk-1)*row+jj, ii) ; 

TE0=temp((kk-1)*row+jj ,ii+1); 

TW0=temp((kk-1)*row+jj,ii-1); 

TN0=temp((kk-1)*row+jj-1, ii); 

TT0=temp((kk-2)*row+jj,ii) ; 

TB0=temp(kk*row+jj, ii) ; 

new(row_num, 2)=apO+teta*(ae+aw+an+2*as*coeffs+at+ab) ; 

new(row_num,3)=-teta*ae; 

new (row__num, 1) =-teta*aw; 

north(row_num,1)=teta*an; 

top (row__num, 1)=teta*at; 

bottom(row_num,1)=teta*ab; 

right(row_num,l)=(l-teta)*(ae*TEO+aw*TWO+an*TNO+at*TTO+ab*TBO)+ ... 
(apO-(1-teta)*(ae+aw+an+2*as*coeffs+at+ab))*TP0 + 2*as*coeffs*Tinf 

2 *as *f luxcoef f s*q_so.uth; 


% for the back surface 

elseif ((kk>l) & (kk<surface) & (jj==l) & (ii>l) & (ii<col)) 
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TPO=temp((kk-1)*row+jj,ii); 

TEO=temp((kk-1)*row+jj,ii+1); 

TWO=temp((kk-1)*row+jj , ii-1); 

TSO=temp((kk-1)*row+jj + l, ii); 

TTO^temp((kk-2)*row+jj, ii); 

TBO=temp(kk*row+j j,ii); 

new(row_num,2)=apO+teta*(ae+aw+2*an*coeffn+as+at+ab); 

new(row^num,3)=-teta*ae; 

new(row_num,1)=-teta*aw; 

south(row_num,1)=teta*as; 

top(row^num,1)=teta*at; 

bottom (row_nuin, 1) =teta*ab; 

right(row_num,l)=(l-teta)*(ae*TE0+aw*TW0+as*TS0+at*TT0+ab*TB0)+ 
(apO-(1-teta)*(ae+aw+2*an*coeffn+as+at+ab))*TP0 + 2*an+coeffn*Tinf 

2*an*f luxcoef fn*q__north; 


% for the left-lateral surface 

elseif ((kk>l) & (kk<surface) & (j j >1) & (jj<row) & (ii==l)) 

TP0=temp((kk-1)*row+jj , ii); 

TE0=temp((kk-1)*row+jj , ii+1); 

TN0=temp((kk-1)*row+jj-1,ii); 

TS0=temp((kk-1)*row+jj + 1, ii) ; 

TT0=temp ( (kk-2) *row+j j , ii) ; 

TB0=temp(kk+row+jj, ii) ; 

new (row__num, 2) =apO+teta* (ae+2 *aw*coef fw+an+as + at + ab) ; 

new (rowjnum, 3) =-teta*ae; 

north(row_num,1)=teta*an; 

south(row_num,1)=teta*as; 

top(row_num,1)=teta*at; 

bottom(row_num,1)=teta*ab; 

right(row_num,1)=(1-teta)*(ae*TE0+an*TN0+as*TS0+at*TT0+ab*TB0)+ ... 
(apO-(1-teta)*(ae+2*aw*coeffw+an+as+at+ab))*TP0 + 2*aw*coeffw*Tinf 

2*aw*fluxcoeffw*q_west; 


% for the right-lateral surface 

elseif ( (kk>l) & (kk<surface) & (j j >1) & (jj<row) & (ii—col) ) 

TP0=temp{(kk-1)*row+jj , ii) ; 

TW0=temp((kk-1)*row+jj,ii-1); 

TN0=temp((kk-1)*row+jj-1, ii); 

TS0=temp((kk-1)*row+j j+l, ii); 

TT0=temp((kk-2)*row+jj,ii); 

TB0=temp(kk*row+j j,ii); 
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new(row_num,2)=apO+teta*(2*ae*coeffe+aw+an+as+at+ab); 

new(row_num,1)=-teta*aw; 

north(row_num,1)=teta*an; 

south(row_num,1)=teta*as; 

top(rowjnum,1)=teta*at; 

bottom(row_num,1)-teta*ab; 

right(row_num,1)=(1-teta)*(aw*TW0+an*TN0+as*TS0+at*TT0+ab*TB0)+ ... 
(apO-(1-teta)*(2*ae*coeffe+aw+an+as+at+ab))*TP0 + 2*ae*coeffe*Tinf 

2*ae*fluxcoeffe*q_east; 


% for the front-top edge 

elseif ( (kk===l) & (jj==row) & (ii>l) & (ii<col)) 

TP0=temp((kk-1)*row+jj,ii); 

TE0=temp((kk-1)*row+jj , ii+1); 

TW0=temp((kk-1)*row+jj,ii-1); 

TN0=temp((kk-1)*row+jj-1,ii); 

TB0=temp(kk*row+jj , ii); 

new(row_num,2)=apO+teta*(ae+aw+an+2*as*coeffs+2*at*coefft+ab); 

new(row_num,3)=-teta*ae; 

new(row_num,1)=-teta*aw; 

north(row_num,1)=teta*an; 

bottom(rowjnum,1)=teta*ab; 

right(row_num,1)=(1-teta)*(ae*TE0+aw*TW0+an*TN0+ab*TB0)+ ... 
(apO-(1-teta)*(ae+aw+an+2*as*coeffs+2*at*coefft+ab))*TP0 + . . . 
(2*as*coeffs+2*at*coefft)*Tinf - 2*as*fluxcoeffs*q_south- 
2*at*fluxcoefft*q_top ... 

-2*at*fluxcoefft*q_up{jj,ii); 
radiation(row^num,1)=-2*at*radcoeff; 
q_boiling_coeff(row^num,1)=-2*at*fluxcoefft; 


% for the front-bottom edge 

elseif ((kk==surface) & (jj==row) & (ii>l) & (ii<col)) 

TP0=temp ( (kk-1) *row+j j , ii) ; 

TE0=temp((kk-1)*row+jj,ii+1); 

TW0=temp((kk-1)*row+jj,ii-1); 

TN0=temp((kk-1)*row+jj-1,ii); 

TT0=temp((kk-2)*row+jj,ii); 

new (row_num, 2) =apO+teta* (ae+aw+‘an+2*as*coef f s+at+2*ab*coef fb) ; 
new(row_num,3)=-teta*ae; 
new(row_num,1)=-teta*aw; 

north(row_num,1)=teta*an; 
top(row_num f 1)=teta*at; 
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right(row_num,1)=(1-teta)*(ae*TEO+aw*TWO+an*TNO+at*TTO)+ ... 
(apO- (1-teta) * (ae+aw+an+2*as*coeffs+at+2*ab*coeffb) ) *TPO + . . 
(2*as*coeffs+2*ab*coeffb)*Tinf - 2*as*fluxcoeffs*q_south- 
2*ab*fluxcoeffb*q_bottom; 


% for the back-top edge 

elseif ( (kk— 1) & (jj==l) & (ii>l) & (ii<col) ) 

TP0=temp((kk-1)*row+j j,ii); 

TE0=temp((kk-1)*row+jj,ii+1); 

TW0=temp((kk-1)*row+jj,ii-1); 

TSO^temp((kk-1)*row+jj + 1,ii); 

TB0=temp(kk*row+j j,ii); 

new(row_num,2)=apO+teta*(ae+aw+2*an*coeffn+as+2*at*coefft+ab); 

new(row_num,3)=-teta*ae; 

new(row^num,1)=-teta*aw; 

south(row_num,1)=teta*as; 

bottom(row_num,l)=teta*ab; 

right(row_num,1) =(1-teta)*(ae*TE0+aw*TW0+as*TS0+ab*TB0)+ . . . 
(apO-(1-teta)*(ae+aw+2*an*coeffn+as+2*at*coefft+ab))*TP0 + .. 
(2*an*coeffn+2*at*coefft)*Tinf - 2*an*fluxcoeffn*q_north- 
2*at*fluxcoefft*q_top... 

-2*at*fluxcoefft*q_up(j j,ii); 
radiation(row_num,1)=-2*at*radcoeff; 
q_boiling_coeff(row_num,1)=-2*at*fluxcoefft; 


% for the back-bottom edge 

elseif ((kk==surface) & (jj==l) & (ii>l) & (ii<col)) 

TP0=temp((kk-1)*row+jj , ii); 

TE0=temp((kk-1)*row+j j,ii+1); 

TWO^temp((kk-1)*row+jj , ii-1) ; 

TS0=temp((kk-1)*row+jj + 1, ii); 

TT0=temp((kk-2)*row+jj , ii); 

new(row_num,2)=apO+teta*(ae+aw+2*an*coeffn+as+at+2*ab*coeffb) ; 

new(row_num,3)=-teta*ae; 

new(row_num,1)=-teta*aw; 

south(row_num,1)=teta*as; 

top(row_num,1)=teta*at; 

right (row_nura, 1) = (1-teta) * (ae*TEO+aw*TWO+as*TSO+at*TTO) + ... 
(apO-(1-teta)* (ae+aw+2*an*coeffn+as+at+2*ab*coeffb))*TP0 + .. 
(2*an*coeffn+2*ab*coeffb)*Tinf - 2*an*fluxcoeffn*q_north- 
2*ab*fluxcoeffb*q_bottom; 
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% for the front-left edge 

elseif ((kk>l) & (kk<surface) & (jj~row) & (ii==l) ) 

TP0=temp ( (kk-1) *row+jj , ii) ; 

TE0=temp ( (kk-1) *row+j j , ii+1) ; 

TN0=temp((kk-1)*row+jj-1,ii); 

TT0=temp((kk-2)*row+jj , ii) ; 

TB0=temp(kk*row+jj,ii); 

new (row__num, 2) =apO+teta* (ae+2*aw*coeffw+an+2*as*coeffs+at+ab) 

new(row_num,3)=-teta*ae; 

north(row_num,1)=teta*an; 

top (row__num, 1) =teta*at; 

bottom(row_num,1)=teta*ab; 

right(row__num,1)=(1-teta)*(ae*TE0+an*TN0+at*TT0+ab*TB0)+ ... 
(apO-(1-teta)*(ae+2*aw*coeffw+an+2*as*coeffs+at+ab))*TP0 + . 
(2*aw*coeffw+2*as*coeffs)*Tinf - 2*aw*fluxcoeffw*q__west- 
2*as*fluxcoeffs*q_south; 


% for the front-right edge 

elseif ((kk>l) & (kk<surface) & (jj==row) & (ii==col)) 

TP0=temp((kk-1)*row+jj,ii); 

TW0=temp((kk-1)*row+jj , ii—1); 

TN0=temp((kk-1)*row+jj-1,ii); 

TT0=temp((kk-2)*row+jj , ii) ; 

TBO-temp(kk*row+jj , ii); 

new(row_num,2)=apO+teta*(2*ae*coeffe+aw+an+2*as*coeffs+at+ab) 

new(row__num, l)=-te'ta*aw; 

north(row_num, 1) ~teta*an; 

top(row^num,1)=teta*at; 

bottom(row_num,1)=teta*ab; 

right(row_num, 1) = (1-teta)*(aw*TW0+an*TN0+at*TT0+ab*TB0)+ ... 
(apO-(1-teta)*(2*ae*coeffe+aw+an+2*as*coeffs+at+ab))*TP0 + . 
(2*ae*coeffe+2*as*coeffs) *Tinf - 2*ae*fluxcoeffe*q_east- 
2*as*f luxcoef fs*q__south; 


% for the back-left edge 

elseif ((kk>l) & (kk<surface) & (jj==l) & (ii==l)) 

TP0=temp((kk-1)*row+jj,ii); 

TE0=temp((kk-1)*row+jj , ii+1); 

TS0=temp((kk-1)*row+jj+1, ii); 

TT0=temp((kk-2)*row+jj,ii); 

TB0=temp(kk*row+j j , ii) ; 
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new(row_num,2)=apO+teta*(ae+2*aw*coeffw+2*an*coeffn+as+at+ab) 

new(row_num,3)=-teta*ae; 

south(row_num,1)=teta*as; 

top(row_num,1)=teta*at; 

bottom(row_num,1)=teta*ab; 

right(row_num,1)=(1-teta)* (ae*TE0+as*TS0+at*TT0+ab*TB0)+ ... 
(apO-(1-teta)* (ae+2*aw*coeffw+2*an*coeffn+as+at+ab))*TP0 + . 
(2*aw*coeffw+2*an*coeffn)*Tinf -2*aw*fluxcoeffw*q_west- 
2*an*fluxcoeffn*q_north; 


% for the back-right edge 

elseif ((kk>l) & (kk<surface) & (jj==l) & (ii==col)) 

TP0=temp((kk-1)*row+jj , ii) ; 

TW0=temp{(kk-1)*row+jj,ii-1); 

TS0=temp((kk-1)*row+jj + 1, ii) ; 

TT0=temp((kk-2)*row+jj , ii) ; 

TB0=temp(kk*row+jj,ii); 

new(row_num, 2) =apO+teta* (2*ae*coef fe+aw+2*an*coef fn+as+at + ab) 

new(row_num,1)=-teta*aw; 

south(row_num,1)=teta*as; 

top'(row_num, 1) -teta*at ; 

bottom(row_num,1)=teta*ab; 

right (row__num, 1) = (1-teta) * (aw*TW0+as*TS0+at*TT0+ab*TB0) + . . . 
(apO-(1-teta)*(2*ae*coeffe+aw+2*an*coeffn+as+at+ab))*TP0 + . 
(2*ae*coeffe+2*an*coeffn)*Tinf - 2*ae*fluxcoeffe*q_east- 
2*an*fluxcoeffn*q_north; 


% for the left-lateral-top edge 

elseif ((kk==l) & (j j >1) & (jj<row) & (ii==l)) 

TP0=temp((kk-1)*row+jj , ii); 

TE0=temp((kk-1)*row+jj,ii+1); 

TN0=temp((kk-1)*row+jj-1, ii) ; 

TS0=temp((kk-1)*row+jj+1,ii) ; 

TB0=temp(kk*row+jj,ii); 

new(row_num,2)=apO+teta*(ae+2*aw*coeffw+an+as+2*at*coefft+ab) 

new(row^num,3)=-teta*ae; 

north(row_num,1)=teta*an; 

south(row^num,l)=teta*as; 

bottom(row_num,1)=teta*ab; 

right(row_num,1) = (1-teta)*(ae*TE0+an*TN0+as*TS0 + ab*TB0)+ . . . 
(apO-(1-teta)*(ae+2*aw*coeffw+an+as+2*at*coefft+ab))*TP0 + . 
(2*aw*coeffw+2*at*coefft)*Tinf - 2*aw*fluxcoeffw*q_west- 
2*at*fluxcoefft*q_top... 

-2*at*fluxcoefft*q_up(jj,ii) ; 
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radiation(row_num,1)=-2*at*radcoeff; 
q_boiling_coeff(row^num,1)=~2*at*fluxcoefft; 


% for the left-lateral-bottom edge 

elseif {(kk==surface) & (jj>1) & (jjcrow) & (ii==l)) 

TP0=temp((kk-1)*row+jj,ii); 

TE0=temp((kk-1)*row+jj,ii+1) ; 

TN0=temp((kk-1)*row+jj-1,ii) ; 

TS0=temp((kk-1)*row+jj+1, ii) ; 

TT0=temp((kk-2)*row+jj , ii) ; 

new(row_num,2)=apO+teta* (ae+2*aw*coeffw+an+as+at+2*ab*coeffb); 

new(row_num,3)=-teta*ae; 

north(row^num,1)=teta*an; 

south(row_num,1)=teta*as; 

top (row__num, 1) =teta*at; 

right(row_num,1)=(1-teta)*(ae*TE0+an*TN0+as*TS0+at*TT0)+ ... 
(apO-(1-teta)*(ae+2*aw*coeffw+an+as+at+2*ab*coeffb))*TP0 + .. 
(2*aw*coeffw+2*ab*coeffb)*Tinf - 2*aw*fluxcoeffw*q_west- 
2*ab* f luxcoef fb*q_bot tom; 


% for the right-lateral-top edge 

elseif ((kk—1) & (j j >1) & (jj<row) & (ii==col)) 

TP0=temp((kk-1)*row+jj , ii) ; 

TWO-temp((kk-1)*row+jj,ii-1) ; 

TN0=temp((kk-1)*row+jj-1, ii) ; 

TS0=temp((kk-1)*row+jj+1, ii) ; 

TB0=temp(kk*row+jj ,ii); 

new(row__num,2)=apO+teta* (2*ae*coeffe+aw+an+as+2*at*coefft+ab); 

new(row_num ; 1)=-teta*aw; 

north(row_num,l)=teta*an; 

south (row__num, 1) =teta*as; 

bottom(row_num,1)=teta*ab; 

right(row_num,1) = (1-teta)* (aw*TWO+an*TNO+as*TSO+ab*TBO)+ ... 
(apO-(1-teta)* (2*ae*coeffe+aw+an+as+2*at*coefft+ab)) *TP0 + .. 
(2*ae*coeffe+2*at*coefft) *Tinf - 2*ae*fluxcoef fe*q__east- 
2*at*fluxcoefft*q__top. . . 

-2*at*fluxcoefft*q_up(jj,ii) ; 
radiation(row_num,l)=-2*at*radcoeff; 
q_boiling__coeff (row_nuin, 1) =-2*at*fluxcoefft; 


% for the right-lateral-bottom edge 

elseif ( (kk—surface) & (jj>l) & (jjcrow) & (ii==col)) 
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TPO=temp((kk-1)*row+jj,ii); 

TWO=temp((kk-1)*row+jj, ii-1) ; 

TNO=temp((kk-1)*row+jj-1,ii) ; 

TSO=temp((kk-1)*row+jj + 1, ii); 

TTO=temp((kk-2)*row+jj, ii) ; 

new(row_num,2)=apO+teta*(2*ae*coeffe+aw+an+as+at+2*ab*coeffb); 

new(row__num, 1) =-teta*aw; 

north(row_num,1)=teta*an; 

south(row_num,1)=teta*as; 

top(row^num,1)=teta*at; 

right(row_num,l)=(l-teta)*(aw*TW0+an*TN0+as*TS0+at*TT0)+ ... 
(apO-(1-teta)*(2*ae*coeffe+aw+an+as+at+2*ab*coeffb))*TP0 + ... 
(2*ae*coeffe+2*ab*coeffb)*Tinf - 2*ae*fluxcoeffe*q_east- 
2*ab*fluxcoeffb*q_bottom; 


% for the interior nodes 

%% as expected, note that the interior node has all of its neighbors 

%% contributing to the coefft matrices 

else 

TP0=temp((kk-1)*row+jj,ii); 

TE0=temp((kk-1)*row+jj,ii+1); 

TW0=temp((kk-1)*row+jj,ii-1) ; 

TN0=temp ( (kk-1) *row+j j-1, ii) ; 

TS0=temp ( (kk-1) *row+j j+1, ii) ; 

TTO-temp ( (kk-2) *row+jj , ii) ; 

TB0=temp(kk*row+jj,ii); 

new(row^num,2)=ap; 
new(row_num,3)=-teta*ae; 
new(row_num,1)=-teta*aw; 
north(row_num,1)=teta*an; 
south(row^num,1)=teta*as; 
top(row_num,1)=teta*at; 
bottom(row^num,1)=teta*ab; 

right(row_num,l)=(l-teta)*(ae*TE0+aw*TW0+an*TN0+as*TS0+at*TT0+ab*TB0)+ 
(apO-(1-teta)* (ae+aw+an+as+at+ab))*TP0; 


end 

% 

the "end" of 

if loop 

to decide where node is located 

end 

% 

the "end" of 

loop for 

ii= . 

. . (cols) 

end 

% 

the "end" of 

loop for 

jj= • 

. . (rows) 

end 

% 

the "end" of 

loop for 

kk= . 

. . (surfaces) 

toe 

% 

end of tic-toc 

loop 
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%% solve the matrix and find the new temperatures by using the 
"aasolvelO .m" 

%% function - row-by-row sweeping iterations from surface to surface 

[temp]=aasolvelO(row,col,surface,new,right,north,south,top,bottom, 

temp, radiation, q_boiling_coef f, Tinf, x_intl, y__intl, x, y) ; 

% print the maximum temp in the entire 3-d grid 
max(max(temp)) 


%% save the row number, column_number, surface number, time and 
temperature 

%% in every iteration 


save_me % print value of savejme counter 

% time_step_interval is the freq at which to save temp data 
time_step_interval=l; 

if (save_me/time_step_interval) =- fix (save__me/time__step_interval) 

%% save__counter is used in "save_thelO.m" function to apend to file¬ 
name 

%% in which temp data is being saved 
save_counter=save_counter+l; 

% save_thelO saves top-surface temp at each time-step into a diff file 
name 

save_thelO (row, col, surface, time, temp (1: row, 1: col) , save__counter, b, x, y) ; 

%% save__a!110 saves all 3-d temp data at each time-step by overwriting 
onto 

%% the same file name - imp for restarting the iterations... 

save_alllO(number,time,temp,b,delta_t,back,front,width,thickness,save_c 
ounter,save_me); 

end 


num_milli_secs=500; % freq at which to save all 3-d temp data 
for aa=l:20 

if (round(savejne) ===round(num_milli_secs*aa) ) 

save_thelO(row,col,surface,time,temp,save_counter,b,x,y); 

end 

end 

end 
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% THE FUNCTION PROGRAM (aasolve 10.m) 

% this program solves the matrix produced by the main program 
% by using iterative sweeping method. 

% 

function [temp] = aasolvelO(row,col,surface,new,right,north,south,top. 


bottom, temp, radiation, q_boiling_coef f, Tinf, x_intl, y__intl, x, y) ; 
% for a a*b*c three dimensional matrix, 


% 


size of new -> (a*b*c)*3 


size of top-> 

size of bottom —> 

size of north -> 

size of south -> 

size of right -> 


a*b*c)*1 
1 
1 
1 
1 


(a*b*c)* 
(a*b*c)* 
(a*b*c)* 
(a*b*c)* 


size of temp.-> (a*b)*c 

size of radiation —> (a*b)’ 
size of q_boiling_coeff —> 
convection heat fluxes from 


TP, TE, TW coefficients for the 
whole grids 

TT coefficients for the whole grids 
TB coefficients for the whole grids 
TN coefficients for the whole grids 
TS coefficients for the whole grids 
The values of the righthandside of 
the matrix 

The temperatures of the grids 
1 The radiation coefficients 
(a*b)*l The boiling and free 
the top surface 


% initialize the temp_old, temp_old is the old temperatures of the 
% grid points. It is used to compare the temperatures before and after 
% the iterations. If the difference is less than 0.1, the iteration 
% is enough. 


temp_old=temp+100*ones(size(temp)); 
it__num=0; 

temp_saturation=100; 

delta__iteration=0.1 ; 
count_iteration=0; 

while max(max(abs(temp_old-temp))) > delta_iteration 
count_iteration = count_iteration+l; 

if count_iteration>deltalteration*50+5 
delta_iteration=delta_iteration+0.1 
end 

max(max(abs(temp_old-temp))) 

it num~it num+1 % iteration number 
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disp (max (max (temp) ) ) 
t emp__o 1 d=t emp; 


% The variables of the boiling regimes 

g=9.81; 
csf-0.0132 ; 
ne=l; 

hfg-2257E3; 
ro_l=957.85; 
ro_ v =0.6; 

surface_tension=58.9E-3; 

Prandtl=l.75; 
miw_l=2.775E-4; 
miw_v=12.02E-6; 
cpl-4211; 
cpv=2029; 
epsi_s=0.82; 
k_liq=0.682; 
k_vap=0.0249; 
alpha_l=0.169E-6; 

T_sat=100; 

T_liquid=27; 
sigma_s=5.67E-8; 


pisi_pisi=zeros(row,col); 


% start the iterations 


for cz=l:surface % for the surfaces- 

for cy=l:row % for the rows 

indexl=(cz-1)*row*col+(cy-1)*col+l; 
index2=(cz-1)*row*col+cy*col; 

newl=new(indexl:index2, : ) ; 
topl=top(indexl:index2,1) ; 
bottoml=bottom(indexl:index2,1); 
northl=north(indexl:index2, 1) ; 
southl=south(indexl:index2,1); 
rightl=right(indexl:index2,1); 
radiationl=zeros(length(indexl:index2),1); 
q_boiling_free=zeros(length(indexl:index2) , 1); 

for cx=l:col % for the columns 


95 



% multiply the coefficients by the temperatures 


if cz==l 

qqqq=0; % heat flux to the outside because of boiling or free 
convection 


% The radiation heat flux from the top surface 

radiationl(l:col,1)^radiation(((cy- 
1)*col+l) : (cy*col),1).*(temp(cy,cx) A 4 ) ; 

radiationl(cx,1)=radiation((cy-1)*col+cx,1)*(temp(cy,cx) A 4); 


% The boiling or free convection heat flux from the top surface 

delta_temp = temp(cy,cx) - temp_saturation; 

if (delta__temp<=5) 

% Rayleigh number 

T_film=(temp(cy,cx)+Tinf)/2+273.15; 
g=9.81; 

beta—0.000000017937*(T_film A 2)+0.0000188*T_film-0.0037515; 
alpha—1.875E-12*(T_film A 2)+1.5525E-9*T_film-1.4995E-7; 
nu=l. 125E-10* (T_film A 2) -8.285E-8*T_f ilm+1 ..5584E-5; 
kconduc=-0.000008125*(T_film A 2)+0.0063775*T_film-0.56895; 

% characteristic length L 

% The free convection coefficient will be found for a 1*1 m2 area 
L=l/4; 


Ra=abs (g*beta* (temp (cy, cx) -Tin*f) * (L A 3) / (nu*alpha) ) ; 
if Ra<16 

Nusselt=l; 


else 


if (Ra<=lE7) 

Nusselt=0.54 *(realpow(Ra,0.25)); 
else 

Nusselt=0.15*(realpow(Ra, (1/3))); 
end 


end 

h_boiling=Nusselt*kconduc/L; 
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q_boiling_free(cx,1)=q_boiling_coeff((cy-1)*col+cx,1)*. . . 

h_boiling*(temp(cy,cx)-Tinf); 

qqqq=h_boiling*(temp(cy,cx)-Tinf); 

end 


if (delta_temp>5 & delta_temp<=30) 

q_prime_s==miw__l*hfg*sqrt (g* (ro_l-ro_v) /surface-tension) *. . . 
(cpl*delta-temp/csf/hfg/Prandtl A ne) A 3; 

tau=pi/3*sqrt(2*pi)*sqrt(surface_tension/g/(ro_l-ro_v))*.. . 
(realpow (abs (ro__v A 2/surface_tension/g/ (ro_l- 

rO-V)),0.25)); 

q^boiling^free (cx, 1) =q_boiling^coeff ( (cy- 
1) *col+cx, 1) *q_prime__s* ... 

(1+(2*k_liq*(T_sat-T—liquid)/. . . 

sqrt(pi*alpha-l*tau))*24/(pi*hfg*ro_v)*... 

realpow(abs(ro^v A 2/surface—tension/g/(ro_l- 

ro_v)),.25)); 

qqqq=q__prime_s*. . . 

(1+(2*k—liq*(T—Sat-T—liquid)/. . . 

sqrt(pi*alpha__l*tau))*24/(pi*hfg*ro_v)*... 

realpow(abs(ro__v A 2/surface_tension/g/(ro_l- 

ro v)),.25)); 


end 

if (delta_temp>30 & delta_temp<=120) 

q_max=0.14 9*hfg*ro_v*realpow(abs (surface_tension*g* (ro__l- 
ro_v)/ro_v A 2),0.25) ; 

q_min=0.09*ro__v*hf g*realpow (abs(surface_tension*g*(ro_l-ro_v)/... 
(ro_l+ro_v) A 2),0.25); 

1 

log_q=logl0 (qjnax/qjmin) /loglO (30/120) *loglO (delta-temp/120)tloglO (q_mi 
n) ; 

q_boiling—free(cx,1)=q_boiling—coeff((cy-1)*col+cx,1)*10 A (log_q) ; 

qqqq=10 A (log_q); 

end 

if (delta_temp>120) 

lambda=2*pi*realpow(abs (surface_tension/g/ (rd—l-ro_v) ), 0.5) ; 
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h_conduct=0.59*realpow(abs(g*(ro_lro_v)*ro_v*k vap A 3*. 
(hfg+O.68*cpv*delta_temp)/(lambda*miw_v*delta_temp)),.25); 

h_radiate=sigma_s*epsi_s*(((temp(cy,cx)) A 4-(T_sat) A 4)/. . . 

(temp (cy, cx) ~T__sat) ) ; 

h_boiling=h_conduct+0.75*h__radiate; 

q_boiling__free (cx, 1) =qJaoiling_coeff ( (cy-1) *col+cx, 1) * . . . 

h_boiling*(temp(cy,cx)-Tinf); 

qqqq=h_boiling*(temp(cy,cx)-Tinf); 
end 

pisi_pisi(cy,cx)=qqqq; 
end 


if cz>l 

topi(cx,1)=topl(cx,1)*temp((cz-2)*row+cy,cx) ; 
end 


if cz<surface 

bottoml(cx,1)=bottoml(cx,1)*temp(cz*row+cy,cx); 
end 

if cy>l 

northl(cx,l)=northl(cx,1)*temp((cz-1)*row+cy-l, cx); 
end 

if cy<row 

southl(cx,1)=southl(cx,1)*temp((cz-1)*row+cy+l, cx); 
end 

end 


% indexl...index2 shows the rows in the new matrix that the 
% calculations are done 


if cz==l 


right1(1:col,1)=rightl(1:col,1)+radiationl(1:col,1)+q_boiling_free(1:co 

1 , 1 ); 

end 

right1(l:col,l)=rightl(l:col,l)+northl(l:col,1) ... 

+southl(1:col,1)+topl(1:col,1)+bottoml(1:col,1); 
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% 


solve the diagonal matrix by using Gauss elimination method 


for i=l:col-1; 

newl(i,3)=newl(i,3)/newl(i, 2); 
right1(i, 1)=rightl(i,1)/newl(i,2); 
newl (i,2)=1; 

newl(i+1,2)-newl(i+1,2)-newl(i+1,1)*newl(i,3) ; 
right1(i+1,1)-right1(i+1,1)-right1(i, 1)*newl(i+1,1); 
newl(i+1,1)=0; 


end 

rightl(col,1)=rightl(col,1)/newl(col,2}; 
newl(col,2)=1; 


for i=col:-l:2 

rightl(i-1,1)=rightl(i-1,1)-rightl(i,1)*newl(i-1,3); 
newl(i-1,3)=0; 

end 


% The calculation is over, update the temperature matrix for the new 
values 

temp((cz-1)*row+cy,:)=rightl(1:col, :) 1 ; 

end 

end 

end 

% figure 

% mesh(x,y,pisi_pisi); 

% pisi_pisi 

% title('heat') 
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% THE FUNCTION PROGRAM (oldjemp 10.m) 

% 

%% The function program (old_templO.m) finds out the new positions of 
%% the moving grid points at the new time step and the old temperature 
%% values of the new grid points by using second degree polynomial 
%% enterpolation. 

function [temp]=old_templO(row,col,surface,temp,x,y,deltax,deltay); 

%% This (for) loop finds out the x-coordinates of the new area where 
%% the grid points move after the time delta_t. 

for forl=l:col 

new_pos_x(1,fori)=x (1,fori)+deltax; 


count = 1; 


while ((x(1,count) < new_pos_x(1,fori)) & (count<col+10)) 
if count==col 

count-col+20; % There was "break" command here before. 

end 

count=count+l; 

end 

if count>=col+20 
count=col; 

end 


if (deltax>0) 

if count>2 % deltax>0 means if the source goes to 

the right 

count — count-2; 
else 

if count > 1 

count =count-l; 

end 

end 

end 

if (deltax<=0) % deltax>0 means if the source goes to 

the right 

if count==col 

count=count-2; 
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else 

if count>l 
count =count-l; 

end 

end 


end 


pos_x(forl) = county- 

end 

%% This (for) loop finds out the y-coordinates of the new area where 
%% the grid points move after the time delta_t. 

for1=0; 

for forl=l:l:row 


new_pos_y(fori)=y(1,fori)tdeltay; 
count = 1; 


while ((y (1, count) > new_pos_y(1,fori)) & (count<row+10)) 

if count-row 
count=row+20; 

end 

count=count+l; 

end 

if count>=row+20 
count=row; 

end 

if (deltay>0) 

if count=-row 
count=count-2; 

else 

if count>l 

count = count~l; 
end 

end 

end 

if (deltay<=0) % deltax>0 means if the source goes to the 

% right 


% deltax>0 means if the source goes to 
% the right 
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if count>2 
count=count-2 ; 
else 

if count==2 
count =count-l; 
end 

end 


end 


po's_y(forl) = count; 


end 

%% Now, we know the coordinates of the area where the grid points 
%% move.We have to use second degree polynomial interpolation to find 
%% the temperature values of the new grid points. Here, we use the 
%% temperature values from the previous time step. 

for for1=1:surface 


for for2=l:row 
for for3=l:col 
% temperatures 

yO=temp((forl-1)*row+pos_y(for2},for3); 
yl=temp((forl-1)*row+pos_y(for2)+1, for3) ; 
y2=temp((forl-1)*row+pos_y(for2)+2, for3) ; 

% positions 

xO=y(pos_y(for2)); 
xl=y(pos_y(for2)+1); 
x2=y(pos_y(for2)+2) ; 

% our second degree polynomial is y=a*x A 2+b*x+c 
% a, b, c are the coefficients. Now find these coefficients... 
b=(yO*(x2 A 2-xl A 2)+yl*(xO A 2-x2 A 2)+y2*(xl A 2-xO A 2) ) / ... 

(xO*(x2 A 2-xl A 2)+xl*(xO A 2-x2 A 2)+x2*(xl A 2-xO A 2) ) ; 

if abs(xO)-=abs(x2) 

a=(y0-y2+b^(x2-x0))/(xO A 2-x2 A 2); 
else 

a= (yO-yl+b*(xl-xO))/(xO A 2-xl A 2); 
end 

c=yO-a*xO A 2-b*xO; 
temp_old((forl- 

1)*row+for2,for3)=a*new_pos_y(for2) A 2+b*new_pos_y(for2)+c; 
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% [pol_coeff]=polyfit(y(pos_y(for2):(pos_y(for2)+2 )), 

% temp((((forl-1)*row+pos_y(for2)):((forl- 

1)*row+pos_y(for2)+2)),for3) *,2); 

% temp_old((forl- 

1) *row+for2, f or3) =polyval (pol_coef f , new_pos__y (for2) ) ; 

end 

end 


end 


for for1=1:surface 


for for2=l:row 
for for3=l:col 


yO=temp_old ( (forl-1) *row+for2, pos__x (for3) ); 
yl=temp_old ( (forl-1) *row+for2, pos_x (for3) +1) ; 
y2=temp_old( (forl-1) *row+for2,pos_x(for3)+2) ; 

xO=x (pos_x (for3) ) ; 

2 $l=x (pos_x (for3) +1) ; 
x2=x (pos__x (f or3)+2) ; 

b=(yO*(x2 A 2-xl A 2)+yl*(xO A 2-x2 A 2)+y2*(xl A 2-xO A 2))/ ... 

(xO*(x2 A 2-xl A 2)+xl*(xO A 2-x2 A 2)+x2*(xl A 2-xO A 2) ) ; 

if abs(xO)~=abs(x2) 

a=(y0-y2+b*(x2-x0))/(xO A 2-x2 A 2) ; 
else 

a== (yO-yl+b* (xl-xO) ) / (xO A 2-xl A 2) ; 
end 

c=yO-a*xO A 2-b*xO; 
temp((forl- 

1) *row+for2, for3) =a*new_pos_x (for3) A 2+b*new_pos_x (for3) +c; 

% [pol_coeff]=polyfit(x(pos_x(for3): (pos_x (for3) + 2 )), ... 

% temp__old( (forl- 

1)*row+for2, ((pos_x(for3)):(pos_x(for3)+2))), 2) ; 

% temp((forl-1)*row+for2,for3)=polyval(pol_coeff,new_pos_x(for3)); 


end 

end 


end 


103 


THE FUNCTION PROGRAM (save alllO.m) 



%% the function program (save_alllO.m) saves all 3-d temp data at each 
%% time-step by overwriting onto 

function 

[]=save_alllO(number,time,temp,b,delta__t, back, front,width,thickness,sav 
e__counter, save_me) ; 

save save alllO 
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% THE FUNCTION PROGRAM (save_thelO.m) 

% 

% the function program (save_thelO) saves top-surface temp at each 
% time-step into a diff file name 

%% (aathlO.m) and (aathlOgoon.m) can not be compiled without putting 
%% "%" symbol in front of saveafter compiling erase "%" symbol and 
%% save the program again. 

function []=save_thelO(row,col,surface,time,temp, save_counter,b,x,y); 
save (['temperaturelO_',num2str(save_counter)]); 
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